I am not sure if it's a mathematics question or coding part although I am pretty sure that the code I used is right,so do excuse me if I a wrong.
I have a set of three first order ode and I am trying to numerically integrate them in python using RK4 method. The problem arises with the exponential term where python rounds it to $0$ and the values thereafter are returned "NAN". I tried various packages to deal with large values still no use. I am thinking if I modify the equations or scale them so that the value is smaller might help, but I am not sure how to do that. Any help or the reference to that would be really helpful as I am kind of in a really stuck up situation. Here's the code I have used for this purpose:
#simple exponential potential # u = K*phi'/H0; v = omega_matter**(1/3)*(1+z); w = l*K*phi' - ln((K**2)*V0/H0**2),z from 1100(initial) to z=0(final) # f,g,h are functions of derivation of u,v,w respectively derived w.r.t t*H0 = T import matplotlib.pyplot as plt import numpy as np import math def f(u,v,w): return -3*u*((v**3 + (u**2)/6 + np.exp(-w)/3)**0.5) + l*np.exp(-w) def g(u,v,w): return -v*(v**3 + (u**2)/6 + np.exp(-w)/3)**0.5 def h(u): return l*u z = np.linspace(start=0.0,stop=1.0,num = 10001) print (z) p = 0.1 q = 1.0 n = 10000.0 dh = (q-p)/n u = [0.0] v =  w = [1.8] l = 1.0 for i in range(0,int(n)): k1 = f(u[i],v[i],w[i]) r1 = g(u[i],v[i],w[i]) s1 = h(u[i]) k2 = f(u[i] + k1*dh/2,v[i] + r1*dh/2,w[i] + s1*dh/2) r2 = g(u[i] + k1*dh/2,v[i] + r1*dh/2,w[i] + s1*dh/2) s2 = h(u[i] + k1*dh/2) k3 = f(u[i] + k2*dh/2,v[i] + r2*dh/2,w[i] + s2*dh/2) r3 = g(u[i] + k2*dh/2,v[i] + r2*dh/2,w[i] + s2*dh/2) s3 = h(u[i] + k2*dh/2) k4 = f(u[i] + dh*k3,v[i] + dh*r3,w[i] + s3*dh) r4 = g(u[i] + k3*dh,v[i] + dh*r3,w[i] + s3*dh) s4 = h(u[i] + dh*k3) u == u.append(u[i] + (dh/6)*(k1 + 2.0*k2 + 2.0*k3 + k4)) v == v.append(v[i] + (dh/6)*(r1 + 2.0*r2 + 2.0*r3 + r4)) w == w.append(w[i] + (dh/6)*(s1 + 2.0*s2 + 2.0*s3 + s4)) plt.plot(z,u, '-b') plt.plot(z,v, '-r') plt.plot(z,w, '-g') #plt.plot(u,v) #plt.plot(u,w) plt.title('quintessence cosmological model') plt.show()