These days, I tried to make a numerical anlysis function for the 4th-order Runge-Kutta method. I could use ODE function in the Scipy or other public libraries, but I wanted to make it and custumize it.

Here is the python code. I refered to someone else's code and textbook.

import numpy as np

def RK4thOrder(func, yinit, x_range, h):
    n = int((x_range[-1] - x_range[0])/h)
    
    x, y = x_range[0], yinit
    
    del x_range
    
    # Containers for solutions
    xsol, ysol = [x,], [y,]

    i = 0
    while i < n:
        k1 = func(x, y)

        yp2 = y + k1 * h / 2

        k2 = func(x+h/2, yp2)

        yp3 = y + k2 * h / 2

        k3 = func(x+h/2, yp3)

        yp4 = y + k3 * h

        k4 = func(x+h, yp4)

        y = y + (((k1 + k4) + (2 * (k2 + k3))) / 6 * h)
        del k1, k2, k3, k4, yp2, yp3, yp4
        x = x + h
        xsol.append(x)
        ysol.append([])
        i = i + 1
        ysol[i] = y
    ysol = np.array(ysol)
    return xsol, ysol

input value:

• func: function of differential equation.
• yinit: inital y value at x[0]. yinit can be 1-dimensional numpy array.
• x_range: [x_min, x_max]
• h: interval for an iteration step.

output value:

• xsol: list of x
• ysol: numpy array of y