Terahertz (THz) spectroscopy is a powerful tool for unambiguously extracting complex-valued material properties (e.g., refractive index, conductivity, etc.) from a wide range of samples, with applications ranging from materials science to biology. However, extracting complex refractive indices from THz time-domain spectroscopy data can prove challenging, especially for multilayer samples. These challenges arise from the large number of transmission-reflection paths the THz pulse can take through the sample layers, leading to unwieldy strings of Fresnel coefficients. This issue has often been addressed using various approximations. However, these approximations are only applicable to specific classes of samples and can give erroneous results when misapplied. An alternative to this approach is to programmatically model all possible paths through the sample. The many paths through the sample layers can be modeled as a tree that branches at every point where the paths diverge, i.e., whenever the pulse can either be transmitted or reflected. This tree can then be used to generate expressions relating the unknown refractive index to the observed time domain data. Here, we provide a freely available open-source package implementing this method as both a MATLAB library and a corresponding graphical user interface, which can also be run without a MATLAB license (https://github.com/YaleTHz/nelly). We have tested this method for a range of samples and compared the results to commonly used approximations to demonstrate its accuracy and wide applicability. Our method consistently gives better agreement than common approximations.