runge-kutta.py

import numpy as np

def runge_kutta(n, y, f, x, h, m, k):
    """
    Perform Runge-Kutta calculation using Gill's method.

    Args:
        n (int): Number of differential equations to be solved.
        y (np.ndarray): Array of N dependent variables with initial values.
        f (np.ndarray): Array of the N derivatives of the variables Y.
        x (float): Independent variable.
        h (float): Step size.
        m (int): Index used in the subroutine.
        k (int): Integer used as the argument of a computed GO TO statement.

    Returns:
        Tuple[int, int]: Updated values of m and k.
    """
    m += 1
    
    if m == 1:
        q = np.zeros(n)
        a = 0.5
    elif m == 2:
        a = 1.707107
        x += 0.5 * h
        q = np.zeros(n)
        for i in range(n):
            y[i] = y[i] + a * (f[i] * h - q[i])
            q[i] = 2.0 * a * h * f[i] + (1.0 - 3.0 * a) * q[i]
        a = 0.2928932
    elif m == 3:
        a = 1.707107
    elif m == 4:
        x += 5 * h
        a = 0.2928932
    
    if a == 0.2928932:
        for i in range(n):
            y[i] = y[i] + h * f[i] / 6.0 - q[i] / 3.0
        m = 0
        k = 2
    
    return m, k

def main():
    """
    Main program to solve differential equations using Runge-Kutta method.

    User provides the initial values for the differential equation, and the
    program calculates the solutions.
    """
    n =  # Set the number of differential equations
    y = np.zeros(n)
    f = np.zeros(n)
    q = np.zeros(n)
    
    x, n_eq, xlim, h, m, *y_values = map(float, input("Enter initial values: ").split())
    y[:n_eq] = y_values
    
    k = 1
    while True:
        if x > xlim:
            break
        
        if k == 1:
            # Programmer must provide expressions for each F() in the main program
            f[0] =  # Add the expression for the 1st derivative
            # Add similar lines for other derivatives if n > 1
            print(f"At X={x}, Y={y[0]}")
        elif k == 2:
            k = 1
            m, k = runge_kutta(n_eq, y, f, x, h, m, k)
        
        x += h

if __name__ == "__main__":
    main()