use the rungekuttagill s method of order 4 to solve ordinary differential equations https www geeksforgeeks orggills4thordermethodtosolvedifferentialequations ravi kumar solve an ordinary differential equations using rungekuttagills method of order 4 args func an ordinary differential equation ode as function of x and y xinitial the initial value of x yinitial the initial value of y stepsize the increment value of x xfinal the final value of x returns solution of y at each nodal point def fx y return xy2 y rungekuttagillsf 0 3 0 2 5 y1 3 4104259225717537 def fx y return x y rungekuttagillsf 1 0 0 2 0 y array 0 0 18 0 32 0 42 0 48 0 5 def fx y return x y y rungekuttagillsf 0 0 0 2 1 traceback most recent call last valueerror the final value of x must be greater than initial value of x def fx y return x y rungekuttagillsf 1 0 0 2 0 traceback most recent call last valueerror step size must be positive solve an ordinary differential equations using runge kutta gills method of order 4 args func an ordinary differential equation ode as function of x and y x_initial the initial value of x y_initial the initial value of y step_size the increment value of x x_final the final value of x returns solution of y at each nodal point def f x y return x y 2 y runge_kutta_gills f 0 3 0 2 5 y 1 3 4104259225717537 def f x y return x y runge_kutta_gills f 1 0 0 2 0 y array 0 0 18 0 32 0 42 0 48 0 5 def f x y return x y y runge_kutta_gills f 0 0 0 2 1 traceback most recent call last valueerror the final value of x must be greater than initial value of x def f x y return x y runge_kutta_gills f 1 0 0 2 0 traceback most recent call last valueerror step size must be positive
from collections.abc import Callable
from math import sqrt

import numpy as np


def runge_kutta_gills(
    func: Callable[[float, float], float],
    x_initial: float,
    y_initial: float,
    step_size: float,
    x_final: float,
) -> np.ndarray:
    
    if x_initial >= x_final:
        raise ValueError(
            "The final value of x must be greater than initial value of x."
        )

    if step_size <= 0:
        raise ValueError("Step size must be positive.")

    n = int((x_final - x_initial) / step_size)
    y = np.zeros(n + 1)
    y[0] = y_initial
    for i in range(n):
        k1 = step_size * func(x_initial, y[i])
        k2 = step_size * func(x_initial + step_size / 2, y[i] + k1 / 2)
        k3 = step_size * func(
            x_initial + step_size / 2,
            y[i] + (-0.5 + 1 / sqrt(2)) * k1 + (1 - 1 / sqrt(2)) * k2,
        )
        k4 = step_size * func(
            x_initial + step_size, y[i] - (1 / sqrt(2)) * k2 + (1 + 1 / sqrt(2)) * k3
        )

        y[i + 1] = y[i] + (k1 + (2 - sqrt(2)) * k2 + (2 + sqrt(2)) * k3 + k4) / 6
        x_initial += step_size
    return y


if __name__ == "__main__":
    import doctest

    doctest.testmod()