Edited version of RK4

RK4-test-edit.c

// edit of this version RK4: https://www.fit.vutbr.cz/study/courses/IMS/public/priklady/rk4-test.c.html

#include <stdio.h>

void Dynamic(double t, double st[] /* *y */, unsigned N, /* out */ double in[] /* *yin */)
{
	// y’’’ + 7y’ - 5 = 0
	in[0] = st[1] * -7 + 5;
	in[1] = st[0];
	in[2] = st[1];
	
//    yin[0] = y[1];
//    yin[1] = -y[0];
}

void RK4step(double t, double state[], unsigned N, double stepsize 
	/* double h, double *y */)
{
    double in[N]; //yin[N];           // integrator input = derivative
    double start[N]; //ystart[N];        // initial state
    // 4 stages results:
    double k1[N];
    double k2[N];
    double k3[N];
    double k4[N];

    int i;
    for (i = 0; i < N; i++)  // save initial value
        start[i] = state[i];

    Dynamic(t, state, in);               // in = f(t, y(t))

    for (i = 0; i < N; i++) {
        k1[i] = h * in[i];     // k1 = h * f(t, y(t))
        state[i] = start[i] + k1[i] / 2;
    }

    Dynamic(t + h / 2, state, in);       // in = f(t+h/2, y(t)+k1/2)

    for (i = 0; i < N; i++) {
        k2[i] = h * in[i];     // k2 = h * f(t+h/2, y(t)+k1/2)
        state[i] = start[i] + k2[i] / 2;
    }

    Dynamic(t + h / 2, state, in);       // in = f(t+h/2, y(t)+k2/2)

    for (i = 0; i < N; i++) {
        k3[i] = h * in[i];     // k3 = h * f(t+h/2, y(t)+k2/2)
        state[i] = start[i] + k3[i];
    }

    Dynamic(t + h, state, in);           // in = f(t+h, y(t)+k3)

    for (i = 0; i < N; i++) {
        k4[i] = h * in[i];     // k4 = h * f(t+h, y(t)+k3)
        // Result:
        // y(t+h) = y(t) + k1/6 + k2/3 + k3/3 + k4/6;
        state[i] = start[i] + k1[i] / 6 + k2[i] / 3 + k3[i] / 3 + k4[i] / 6;
    }
    // Return: y = new state at time t+h
}

// Experiment control
int main()
{
	const int N = 3; 
    double y[N] = { 0, 0, 1 };
    double stepsize = 0.1;

    double t = 0;
    while (t < 20) {
        printf("%8.3f % -11g % -11g % -11g  % e\n", t, y[0], y[1], y[2], y[0] - sin(t));
        RK4step(t, y, N, stepsize);        // make step
        t += stepsize;          // advance the simulation time
    }
    printf("%8.3f % -11g % -11g % -11g  % e\n", t, y[0], y[1], y[2], y[0] - sin(t));
    return 0;
}