# Solving perturbed Einstein Boltzmann equations using RK4

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$$

• If you can rearrange everything into a system of the form $\dot{y} = f(y)$, then you can just apply the method from your book or Wikipedia or something. It seems like most of the work is to properly substiture all the time derivatives to get an ODE for each component. Commented Aug 24, 2023 at 15:12
• Why don't you start with a simpler case -- say, a single ODE like $\dot x=3x$ -- and then move on to a system of two simple ODEs. If you don't even know where to start, then this isn't your time to look at such a complex system yet. Commented Aug 24, 2023 at 17:06