The quantitative integration of nuclear measurements into the in situ petrophysical and geophysical evaluation of rock formations has been elusive because of the lack of efficient algorithms to simulate them. We discovered a new method for rapid numerical simulation of borehole neutron measurements using Monte Carlo (MC)-derived spatial flux sensitivity functions (FSFs) and diffusion flux-difference (DFD) approximations. The method calculates spatial sensitivity flux perturbations using flux-difference approximations of one-group neutron diffusion models. By invoking appropriate boundary conditions, the one-group, time-independent neutron diffusion solution is implemented for nonmultiplying systems in 2D and 3D cylindrical coordinates. The solution is differentiated with respect to the neutron cross section to obtain an expression for flux-difference due to cross-section perturbations. Constant transport-correction coefficients for cross-section parameters are calculated with a flux-fitting method to account for deviations of borehole neutron measurements from the physics of diffusion. Thereafter, spatial FSF responses are rapidly and accurately calculated using a first-order Rytov DFD approximation. Estimated flux-differences are next used to calculate lumped higher order perturbation terms. The DFD technique is tested and benchmarked against MC calculations in the presence of standoff, invasion, and well deviation for wireline and logging-while-drilling tools. Benchmark examples and application in highly deviated wells indicate that neutron porosity logs can be accurately and efficiently simulated with the new DFD method, even in complex geometrical and physical conditions, with errors lower than 1.2 porosity units.