Stellar Nucleosynthesis

I've been working the last few months in Stellar Nucleosynthesis, more specific in the proton-proton chain decay, which is considered the most important nuclear process given in the Sun because is the beginning of a series of nuclear reactions and also it takes the longest period of the star's life. Basically I used the four principal nuclear reactions that occur in the Sun where the product of one reaction is the starting material of the next reaction, leading from Hydrogen to Helium. Using some approximations in order to obtain a simpler model we got only three reactions.From these reactions we can get the differential equations that describe the variation of the number of nuclei with time. Using the basic model for the velocities of the reactions, in which we are no taking account the temperature, just the concentration of nuclei, we get this system of ODE's \begin{eqnarray} \frac{d[H]}{dt}&=&k_{-1}[D]-k_1[H]^2+k_{-2}[^3He]-k_2[D][H]\nonumber\\ &+&k_3[^3He]^2-k_{-3}[^4He][H]^2\\ \frac{d[D]}{dt}&=&k_1[H]^2-k_{-1}[D]+k_{-2}[^3He]-k_2[D][H]\\ \frac{d[^3He]}{dt}&=&k_2[D][H]-k_{-2}[^3He]+k_{-3}[^4He][H]^2-k_3[^3He]^2\\ \frac{d[^4He]}{dt}&=&k_3[^3He]^2-k_{-3}[^4He][H]^2 \end{eqnarray} That we solved it doing these ODE's dimensionless and using the numerical Euler's method.