Program for RK-4 first order IVP

Classical RK-4 for dy/dx = f(x, y), y(x₀) = y₀, step h:

k₁ = h·f(xₙ, yₙ)
k₂ = h·f(xₙ + h/2, yₙ + k₁/2)
k₃ = h·f(xₙ + h/2, yₙ + k₂/2)
k₄ = h·f(xₙ + h, yₙ + k₃)
yₙ₊₁ = yₙ + (k₁ + 2k₂ + 2k₃ + k₄)/6

Algorithm

1. Read x₀, y₀, the end point xₙ and the step h.
2. Compute the number of steps n = (xₙ − x₀)/h.
3. Repeat n times: find k₁, k₂, k₃, k₄, update y and x, print them.

C program

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

/* dy/dx = f(x, y): change this function for another problem */
double f(double x, double y)
{
    return x + y;
}

int main(void)
{
    double x0, y0, xn, h, x, y, k1, k2, k3, k4;
    int i, n;

    printf("Enter x0, y0, xn and step size h: ");
    scanf("%lf %lf %lf %lf", &x0, &y0, &xn, &h);

    n = (int)((xn - x0) / h + 0.5);      /* number of steps */
    x = x0;
    y = y0;

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

    for (i = 0; i < n; i++)
    {
        k1 = h * f(x, y);
        k2 = h * f(x + h / 2, y + k1 / 2);
        k3 = h * f(x + h / 2, y + k2 / 2);
        k4 = h * f(x + h, y + k3);

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

Sample run

For f = x + y, input 0 0 0.4 0.1 (y(0) = 0, up to x = 0.4, h = 0.1):

  x         y
  0.0000      0.0000
  0.1000      0.0052
  0.2000      0.0214
  0.3000      0.0499
  0.4000      0.0918

The exact value y = eˣ − x − 1 gives y(0.4) = 0.0918. To solve another problem, change only the body of f.