# ODE45 solver implementation in Python

I have a project where I need ODE solver without dependencies to libraries like Scipy. I decide to implement ODE45. According to tutorials from internet and from what I remember from classes I implement it somehow.

The code below contains also example function (pendulum) with exactly same values as in the scipy ODE tutorial https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.integrate.odeint.html. The resulting plot from my code looks the same.

However, I am not really sure if I done it in the good way. Especially I would really appreciate any hints how to improve the functions ode45_step and ode45 and verification whether my implementation is correct. I do not need to check the input variables, it will be done somewhere else.

import numpy as np
import matplotlib.pylab as plt

def ode45_step(f, x, t, dt, *args):
"""
One step of 4th Order Runge-Kutta method
"""
k = dt
k1 = k * f(t, x, *args)
k2 = k * f(t + 0.5*k, x + 0.5*k1, *args)
k3 = k * f(t + 0.5*k, x + 0.5*k2, *args)
k4 = k * f(t + dt, x + k3, *args)
return x + 1/6. * (k1 + k2 + k3 + k4)

def ode45(f, t, x0, *args):
"""
4th Order Runge-Kutta method
"""
n = len(t)
x = np.zeros((n, len(x0)))
x[0] = x0
for i in range(n-1):
dt = t[i+1] - t[i]
x[i+1] = ode45_step(f, x[i], t[i], dt, *args)
return x

def f(t, y, b, c):
"""
Pendulum example function.
"""
theta = y[0]
omega = y[1]
dydt = [omega, -b*omega - c*np.sin(theta)]
return np.array(dydt)

b = 0.25
c = 5.0

N = 101

x0 = np.array([np.pi - 0.1, 0.0])
t = np.linspace(0, 10, N)

x = ode45(f, t, x0, b, c)

plt.plot(x)
plt.show()


return x + 1/6. * (k1 + k2 + k3 + k4)

return x + 1/6. * (k1 + 2*k2 + 2*k3 + k4)