#include <stdio.h>
#include <stdlib.h>
void rk(unsigned n, void fun(unsigned, double, double*, double*),
double h, double x, double* y, double* f, double* s, double* yk)
{
int i;
double xk;
//double f[N_EQ], s[N_EQ], yk[N_EQ];
fun(n, x, y, f);
for(i = 0; i < n; ++i)
{ s[i] = f[i]; yk[i] = y[i] + (h/2.)*f[i]; }
xk = x + h/2.;
fun(n, xk, yk, f);
for(i = 0; i < n; ++i)
{ s[i] = s[i] + 2.*f[i]; yk[i] = y[i] + (h/2.)*f[i]; }
fun(n, xk, yk, f);
for(i = 0; i < n; ++i)
{ s[i] = s[i] + 2.*f[i]; yk[i] = y[i] + h*f[i]; }
xk = x + h;
fun(n, xk, yk, f);
for(i = 0; i < n; ++i)
{ y[i] = y[i] + (h/6.)*(s[i] + f[i]); }
}
