/* fitz.c  FitzHugh-Nagumo equations  system of two ODE's */

#include <stdio.h>
#undef  abs
#define abs(x) ((x)<0.0?(-(x)):(x))

double vpf(double v, double r) /* independent variable t not needed */
{
  return v - v*v*v/3.0 + r;
}

double rpf(double v, double r) /* independent variable t not needed */
{
  return -(v - 0.2 - 0.2*r);
}


int main(int argc, char argv[])
{
  double v, r, vp, rp, h, t, tn, step;
  double v1, v2, vp2, r1, r2, rp2;
  int i;

  printf("fitz.c running, solve: \n");
  printf(" v' = v - v^3/3.0 + r      initial v = -1.0 \n");
  printf(" r' = -(v - 0.2 - 0.2*r)   initiaL R = 1.0 \n");

  h = 0.001;
  t = 0.0;
  step = 0.25;
  tn = step;
  v = -1.0;
  r = 1.0;
  while(t<25.0)
  {
    vp = vpf(v,r);
    rp = rpf(v,r);
    v1 = v + h * vp;
    r1 = r + h * rp;
    v2 = v + h/2.0 * vp;
    r2 = r + h/2.0 * rp;
    vp2 = vpf(v2, r2);
    rp2 = rpf(v2, r2);
    v2 = v2 + h/2.0 * vp2;
    r2 = r2 + h/2.0 * rp2;
    if(t==0.0)
    {
      printf("t=%g, v=%g, v1=%g, v2=%g, vp=%g \n", t, v, v1, v2, vp);
      printf("t=%g, r=%g, r1=%g, r2=%g, rp=%g \n", t, r, r1, r2, rp);
    }
    if(abs(v1-v2) + abs(r1-r2) > 0.000001)
    {
      h = h/2.0;
      printf("h reduced to %g \n", h);
      if(h<1.0e-6) break;
      continue;
    }
    v = v2;
    r = r2;
    t = t + h;
    if(t>=tn)
    {
      printf("t=%g, vp=%g, rp=%g, v=%g, r=%g \n", t, vp, rp, v, r);
      tn = tn + step;
    }   
  }

  return 0;
} /* end fitz.c */