We present a method which allows for the extraction of physical quantities directly from zero- to ultralow-field nuclear magnetic resonance (ZULF NMR) data. A numerical density matrix evolution is used to simulate ZULF NMR spectra of several molecules in order to fit experimental data. The method is utilized to determine the indirect spin-spin couplings ($J$-couplings) in these, which is achieved with precision of $10^{-2}$--$10^{-4}$ Hz. The simulated and measured spectra are compared to earlier research. Agreement and precision improvement for most of the $J$-coupling estimates are achieved. The availability of an efficient, flexible fitting method for ZULF NMR enables a new generation of precision-measurement experiments for spin-dependent interactions and physics beyond the Standard Model. ; Comment: 8 pages, 7 figures