I'm trying to learn to numerically solve the perturbed Boltzmann-Einstein equations in cosmology using the RK4 method. These are the equations:
$$\dot{\Theta}_{r,0}+k\Theta_{r,1}=-\dot{\Phi}$$
$$\dot{\Theta}_{r,1}+\frac{k}{3}\Theta_{r,0}=\frac{-k}{3}\Phi $$
$$\dot{\delta}+ikv=-3\dot{\Phi}$$
$$\dot{v}+\frac{\dot{a}}{a}v=ik\Phi $$
$$\dot{\Phi}=\frac{1}{3\dot{a}}\frac{3H_{0}^{2}}{2}(\Omega_{m}\delta+4\Omega_{r}\Theta_{r,0}a^{-1})-ak^{2}\Phi-\frac{\dot{a}}{a}\Phi $$
How do I go about this?