#include "f.hh"

void rk4(double a, double b, double h, int N, double y0, double rkout[])
{
  double y=y0;
  double x=a;
  double k1,k2,k3,k4;

  for(int i=0; i<N; i++)
    {
      rkout[i]=y;
      k1=h*f(x,y);
      k2=h*f(x+(h/2.0), y+(k1/2.0));
      k3=h*f(x+(h/2.0), y+(k2/2.0));
      k4=h*f(x+h, y+k3);
      y = y + (k1/6.0) + (k2/3.0) + (k3/3.0) + (k4/6.0);
      x=x+h;
    }
}