#include "f.hh"
#include "math.h"

//
// Heun function
//
void heun (double a, double b, double h, int N,double y0, double heunout[])
{
  double y=y0;
  double x=a;
  double k1, k2;
  double temp;

  for(int i=0; i<N; i++)
    {
      heunout[i] =y;
      k1 =f(x,y);
      temp = y + h * k1;
      k2 = f(x+h, temp); 
      y = y + (h/2.) * (k1 + k2);
      x=x+h;
    }
}