Numerically solving scalar field model of dark energy Not sure if I should be asking this in maths or physics section so please excuse me if asked in the wrong section. I am trying to numerically integrate the equations obtained in the exponential potential scalar field model of dark energy. The set of equations I have obtained is:
\begin{align}
x &= tH_{0} \\
y_{1} &= \frac{a}{a_{0}} \\
y_{2} &= \frac{\psi}{\psi _{0}}\\
y_{3} &= \frac{dy_{2}}{dx}
\end{align}
which on differentiating w.r.t x I get:
\begin{align}
{y_{1}}\prime &= y_{1}\left(\frac{\Omega_{m}}{y_{1}^3} + \frac{y_{3}^2}{6.0} + \frac{V(y_{2})}{2.H_{0}^2}\right) \\
y_{2}\prime &= y_{3} \\
y_{3}\prime &= -3\frac{y_{1}\prime}{y_{1}}y_{3} - \frac{1}{H_{0}^2}\frac{\partial v(y2)}{\partial y_{2}}
\end{align}
My doubt here is I am using RK4 method for solving this and when I write the increment should I increment the independent variable $x$ too?Like when I write formulae for $k_1,\,k_2,\,k_3,\,k_4$ etc I am incrementing $y_1,\,y_2,\,y_3$. As far as I know it is not needed as these ODE are autonomous equations yet when I try to plot them I am not getting the expected result.  
Another doubt is I am plotting it w.r.t scale factor( $a=0$ to $a=1$) and when I plot $y_1$ which is $y_{1} = \frac{a}{a_{0}}$ shouldn't I be getting a linear curve or am I missing on some major concept? 

 A: Let's look at what $y_1 = \frac{a}{a_0}$ means.
Since $a = a_0 e^N$, where $N$ is the e-fold number, $\frac{a}{a_0} = e^N$ and so $y_1 = e^N$.
When you plot $y_1$ against $N$, you do indeed get the green curve shown in the plot in your question (using a calculator to plug in values of $N$ to get $e^N$ gives a rough approximation).
I think the next step would be for you to plot the Hubble parameter $H$ against redshift $z$ (you can find $H$ by using $H = \frac{\dot{a}}{a}$). Depending on the values you are using for the slope of the potential (generally called $\lambda$ in the literature) and the initial values of the potential $V_0$ and the field itself ($\phi$ and $\dot{\phi}$)*, you should be able to find a value of $H_0$ consistent with observational data 
i.e. $H_0 \approx 70$ kms$^{-1}$Mpc$^{-1}$.
Put $\Omega_m$ and $\Omega_\phi$ on a separate plot too. If you plot these against redshift you should again be able to find a set of model parameters where $\Omega_m \approx 0.3$ and $\Omega_\phi \approx 0.7$, again in line with current observations.
Also, if you want to make life a little easier for yourself, use from scipy import integrate and then integrate.odeint. This is a nice function in Python and saves you from having to solve the equations numerically.

*I was coupling the field to matter, but found that values of
$\lambda = 0.23$, $V_0 = 0.76 \times 10 ^{-120}$, $\phi_0= 0.7$ and $\dot{\phi_0}= 0$ worked nicely. I was matching the Planck data, so this gave $H_0 = 67.8$ kms$^{-1}$Mpc$^{-1}$.
