Program for RK-4 second order IVP

A second-order IVP y'' = g(x, y, y') with y(x₀) = y₀, y'(x₀) = z₀ is reduced to two first-order equations by putting z = y':

dy/dx = z            = f₁(x, y, z)
dz/dx = g(x, y, z)   = f₂(x, y, z)

RK-4 is then applied to both equations together:

k₁ = h·f₁(x, y, z)                     l₁ = h·f₂(x, y, z)
k₂ = h·f₁(x+h/2, y+k₁/2, z+l₁/2)       l₂ = h·f₂(x+h/2, y+k₁/2, z+l₁/2)
k₃ = h·f₁(x+h/2, y+k₂/2, z+l₂/2)       l₃ = h·f₂(x+h/2, y+k₂/2, z+l₂/2)
k₄ = h·f₁(x+h, y+k₃, z+l₃)             l₄ = h·f₂(x+h, y+k₃, z+l₃)
y ← y + (k₁+2k₂+2k₃+k₄)/6              z ← z + (l₁+2l₂+2l₃+l₄)/6

C program (example: y'' = y − x·y')

#include <stdio.h>
#include <math.h>

/* y'' = g(x, y, z) with z = y'.  Example: y'' = y - x*y' */
double f1(double x, double y, double z)      /* dy/dx = z */
{
    return z;
}
double f2(double x, double y, double z)      /* dz/dx = y'' */
{
    return y - x * z;
}

int main(void)
{
    double x0, y0, z0, xn, h, x, y, z;
    double k1, k2, k3, k4, l1, l2, l3, l4;
    int i, n;

    printf("Enter x0, y(x0), y'(x0), xn and h: ");
    scanf("%lf %lf %lf %lf %lf", &x0, &y0, &z0, &xn, &h);

    n = (int)((xn - x0) / h + 0.5);
    x = x0; y = y0; z = z0;

    printf("\n  x         y          y'\n");
    printf("%8.4f  %10.4f  %10.4f\n", x, y, z);

    for (i = 0; i < n; i++)
    {
        k1 = h * f1(x, y, z);
        l1 = h * f2(x, y, z);
        k2 = h * f1(x + h / 2, y + k1 / 2, z + l1 / 2);
        l2 = h * f2(x + h / 2, y + k1 / 2, z + l1 / 2);
        k3 = h * f1(x + h / 2, y + k2 / 2, z + l2 / 2);
        l3 = h * f2(x + h / 2, y + k2 / 2, z + l2 / 2);
        k4 = h * f1(x + h, y + k3, z + l3);
        l4 = h * f2(x + h, y + k3, z + l3);

        y = y + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
        z = z + (l1 + 2 * l2 + 2 * l3 + l4) / 6;
        x = x + h;
        printf("%8.4f  %10.4f  %10.4f\n", x, y, z);
    }
    return 0;
}

Sample run

Input 0 1 0 0.2 0.1 (y(0) = 1, y'(0) = 0, up to x = 0.2, h = 0.1):

  x         y          y'
  0.0000      1.0000      0.0000
  0.1000      1.0050      0.0998
  0.2000      1.0199      0.1987

For any other equation change only the f2 function (the right side y'' = …).