#include "stuff.h"

/* Given values for variables y[0..M-1] at x, this function implements */
/* a fourth order RK scheme to advance them to yout[0..M-1] over */
/* a step size dx. */
/* Need to supply a function f which computes the RHS of the
/* equations dy_i/dx=f_i(x,y) and returns the dy[M] vector i=0..M-1

#define M 2

void rk4(double y[], double yout[], double x, double dx){
double k1[M],k2[M],k3[M],k4[M];
double yt[M],dy[M];
int i;

f(x,y,dy);
for(i=0;i<M;i++){
k1[i]=dx*dy[i];
yt[i]=y[i]+0.5*k1[i];}

f(x+dx*0.5,yt,dy);
for(i=0;i<M;i++){
k2[i]=dx*dy[i];
yt[i]=y[i]+0.5*k2[i];}

f(x+dx*0.5,yt,dy);
for(i=0;i<M;i++){
k3[i]=dx*dy[i];
yt[i]=y[i]+k3[i];}

f(x+dx,yt,dy);
for(i=0;i<M;i++)
k4[i]=dx*dy[i];

for(i=0;i<M;i++)
yout[i]=y[i]+(1.0/6.0)*(k1[i]+2*k2[i]+2*k3[i]+k4[i]);

return;
}