from numba import jit
from scipy.integrate import quad
import numpy as np

@jit(nopython = True)
def integrand1(x,e,delta,r):
    return (-2*np.sqrt(e*r)/np.pi)*(x/np.sqrt(1-x**2))/(1+(delta+2*x*np.sqrt(e*r))**2)

@jit(nopython = True)
def f1(e,delta,r):
    return quad(integrand1, -1, 1, args=(e,delta,r))[0]

@jit(nopython = True)
def runge1(e,dtau,delta,r):
    k1 = f1(e,delta,r)
    k2 = f1((e+k1*dtau/2),delta,r)
    k3 = f1((e+k2*dtau/2),delta,r)
    k4 = f1((e+k3*dtau),delta,r)
    return e + (dtau/6)*(k1+2*k2+2*k3+k4)

time_steps = 60
e = 10
dtau=1
r=1
delta=-1

e_array = np.zeros(time_steps)
time = np.zeros(time_steps)
for i in range(time_steps):
    e_array[i] = e
    time[i] = i*dtau
    e = runge1(e,dtau,delta,r)

e_array = pythoncode.py(time_steps = 60, e = 10, dtau = 1, r = 1, delta = -1)