exponential-rk4

import matplotlib.pyplot as plt
import numpy as np
import math


def f(y1,y2,y3):
    return y1*(omega_matter/(y1**3) + (H**2)*(y3**2)/(2.0*rho_c) + (v0*np.exp(-l*y2*k))/(rho_c))**(1.0/2.0)

def g(y3):
    return y3

def h(y1,y2,y3):
    return -3*y3*(H**2)*(omega_matter/(y1**3) + (H**2)*(y3**2)/(2.0*rho_c) + (v0*np.exp(-l*y2*k))/(rho_c))**0.5 + (l*k*v0*np.exp(-l*y2*k))/(H**2)


a = np.linspace(start=1.0,stop=10**(-3),num=10001)
N = np.log(a)

#constants used
H = 2.27e-18
l = 1.5
G = 6.637*(10**(-11))
k = (8*3.14*G)**0.5
om_de = 0.75
omega_matter = 1 - om_de
w0 = -0.85
rho_c = 3*(H**2)/(k**2)
c = ((om_de*(1 + w0)*rho_c)/(H**2))**0.5
v0 = (om_de*(1 - w0)*rho_c)/(2*np.exp(-l*k))

#initial conditions
y1 = [1.0]
y2 = [1.0]
y3 = [c]

#the parameters for integration
p = 0.0
q = 1.0
n = 10000.0
dh = -(q-p)/n

for i in range(0,int(n)):
    k1 = f(y1[i],y2[i],y3[i])
    r1 = g(y3[i])
    s1 = h(y1[i],y2[i],y3[i])
    k2 = f(y1[i] + k1*dh/2.0,y2[i] + r1*dh/2.0,y3[i] + s1*dh/2.0)
    r2 = g(y3[i] + s1*dh/2.0)
    s2 = h(y1[i] + k1*dh/2.0,y2[i] + r1*dh/2.0,y3[i] + s1*dh/2.0)
    k3 = f(y1[i] + k2*dh/2.0,y2[i] + r2*dh/2.0,y3[i] + s2*dh/2.0)
    r3 = g(y3[i] + s2*dh/2.0)
    s3 = h(y1[i] + k2*dh/2.0,y2[i] + r2*dh/2.0,y3[i] + s2*dh/2.0)
    k4 = f(y1[i] + dh*k3,y2[i] + dh*r3,y3[i] + dh*s3)
    r4 = g(y3[i] + dh*s3)
    s4 = h(y1[i] + dh*k3,y2[i] + dh*r3,y3[i] + dh*s3)
    y1 == y1.append(y1[i] + (dh/6.0)*(k1 + 2.0*k2 + 2.0*k3 + k4))
    y2 == y2.append(y2[i] + (dh/6.0)*(r1 + 2.0*r2 + 2.0*r3 + r4))
    y3 == y3.append(y3[i] + (dh/6.0)*(s1 + 2.0*s2 + 2.0*s3 + s4))

#other calculations
A = [(H)*((omega_matter/(y1[i]**3) + (H**2)*(y3[i]**2)/(2.0*rho_c) + (v0*np.exp(-l*y2[i]*k))/(rho_c))**(1.0/2.0)) for i in range(len(y3))]
B = [y2[i]/(10**4) for i in range(len(y2))]
C = [y3[i]/(10**4) for i in range(len(y3))]

w = [((H**2)*y3[i]**2 - 2.0*v0*np.exp(-l*k*y2[i]))/((H**2)*y3[i]**2 + 2.0*v0*np.exp(-l*k*y2[i])) for i in range(len(y3))]
omde = [(((H**2)*(y3[i]**2))/2.0 + (v0*np.exp(-l*k*y2[i]))/rho_c) for i in range(len(y3))]
omnr = [(1 - omde[i]) for i in range(len(omde))]
phi_dot = [y3[i]*H for i in range(len(y2))]

#plt.plot(N,phi_dot,'y',label='phi_dot')
#plt.plot(N,omnr,'-r',label='matter density')
#plt.plot(N,omde,'-b',label='scalar density')
plt.plot(N,w,'-k',label='eos')
#plt.plot(N,B,'-g',label='y2')
#plt.plot(N,C,'-y',label='y3')
#plt.plot(a,A,'-y',label='hubble parameter')
#plt.plot(N,y1,'-g',label='y1')
#plt.plot(N,y2,'-m',label='y2')
#plt.plot(N,y3,'-y',label='y3')
plt.xlabel('e-foldings')
plt.title('exponential potential-b')
plt.legend(loc='best')
plt.show()