A response matrix method (RMM) is applied to Monte Carlo simulations to efficiently compute neutron pulse height distributions (PHDs) in organic scintillation detectors. The PHD calculations and their associated uncertainty are compared for a polyethylene-shielded and lead-shielded 252Cf source for three different techniques: fully analog MCNPX-PoliMi, the RMM, and the RMM with source biasing. The RMM with source biasing reduces computation time or improves the figure of merit on average by a factor of 600 for polyethylene shielding and a factor of 300 for lead shielding (when compared to the fully analog calculation). The simulated neutron PHDs show good agreement with the laboratory measurements, thereby validating the RMM.