We present the use of the Douglas-Gunn Alternating Direction Implicit finite difference method for computationally efficient simulation of the electric field propagation through a wide variety of optical fiber geometries. The method can accommodate refractive index profiles of arbitrary shape and is implemented in a tool called BPM-Matlab. We validate BPM-Matlab by comparing it to published experimental, numerical, and theoretical data and to commercially available state-of-the-art software. It is user-friendly, fast, and is available open-source. BPM-Matlab has a broad scope of applications in modeling a variety of optical fibers for diverse fields such as imaging, communication, material processing, and remote sensing.