Numerical simulations of neutron production from deuterium–lithium nuclear fusion reactions have been performed. A set of differential cross sections for the 7Li(d,xn) reaction for incident deuteron energies of up to 50 MeV is assembled. The angular distribution of neutrons from a thick lithium target is simulated and benchmarked against experimental data. Two-stage neutron production from laser–target experiments has been studied as a function of laser intensity and energy. During the first stage a well collimated deuteron beam is generated using a high-intensity ultrashort pulse laser. During the second stage it is transported through a lithium target using a 3D Monte-Carlo ion beam–target deposition model. The neutron yield is estimated to be ∼108 neutrons J−1 laser energy. Some 1010 neutrons can be expected from a ∼100 J petawatt-class laser. For incident deuteron energies above 1 MeV the proposed scheme for neutron production from d–Li reactions is superior to that from d–d reactions, producing a collimated beam of neutrons with higher neutron yield.