The first three moments of the energy distributions of products from fusion reactions in thermonuclear plasmas with Maxwellian ion velocity distributions are determined analytically. Relativistic kinematics is used allowing the desired accuracy to be reached, which is 2 to 3 orders of magnitude better than previous analytical results. In particular, neutron spectra of the reactions D(d,n)3He and D(t,n)α for plasma ion temperatures 0 < Ti < 100 keV are studied for which the results are also given in tabulated form with interpolation formulas for the purpose of practical use. The neutron energy distributions are also calculated with numerical methods, in order to assess the analytical results. High accuracy calculations are motivated by the crucial role that neutron measurements are envisaged to play in the next step fusion experiments on burning plasmas, which are also discussed.