We present a novel technique to numerically solve pulsed optical beam or “light bullet” propagation in bulk nonlinear media based on the scalar nonlinear Schrödinger (NLS) equation in three dimensions with spherical symmetry. Using fast algorithms for spherical Fourier Bessel transforms along with adaptive longitudinal stepping and transverse grid management in a symmetrized split-step technique, it is possible to accurately study many nonlinear effects, including the possibility of spatio-temporal collapse, or the collapse-arresting mechanism due to saturable nonlinearity.