/* fitz2.c  FitzHugh-Nagumo equations  system of two ODE's     */
/*          using Runge-Kutta-Fehlberg method similar to ode45 */

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

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

double rpf(double t, double v, double r) /* independent variable t not used */
{
  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;
  double k1, k2, k3, k4, k5, k6;
  int i;

  printf("fitz2.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)
  {
    k1 = h * vpf(t,v,r); /* full implementation, even though some not used */
    k2 = h * vpf(t+h/4.0,v+k1/4.0,r);
    k3 = h * vpf(t+3.0*h/8.0,v+3.0*k1/32.0+9.0*k2/32.0,r);
    k4 = h * vpf(t+12.0*h/13.0,v+1932.0*k1/2197.0-7200.0*k2/2197.0+
                 7296*k3/2197.0,r);
    k5 = h * vpf(t+h,v+439.0*k1/216.0-8.0*k2+3680.0*k3/513.0-
                 845.0*k4/4104.0,r);
    k6 = h * vpf(t+h/2.0,v-8.0*k1/27.0+2.0*k2-3544.0*k3/2565.0+
                 1859.0*k4/4104.0-11.0*k5/40.0,r);
    v1 = v + (25.0*k1/216.0+1408.0*k3/2565.0+2197.0*k4/4104.0-k5/5.0);
    vp =     (16.0*k1/135.0+6656.0*k3/12825.0+28561.0*k4/56430.0-
              9.0*k5/50.0+2.0*k6/55.0);
    vp = vp/h;    /* for printout */
    v2 = v + h * vp;

    k1 = h * rpf(t,v,r);
    k2 = h * rpf(t+h/4.0,v,r+k1/4.0);
    k3 = h * rpf(t+3.0*h/8.0,v,r+3.0*k1/32.0+9.0*k2/32.0);
    k4 = h * rpf(t+12.0*h/13.0,v,r+1932.0*k1/2197.0-7200.0*k2/2197.0+
                 7296*k3/2197.0);
    k5 = h * rpf(t+h,v,r+439.0*k1/216.0-8.0*k2+3680.0*k3/513.0-
                 845.0*k4/4104.0);
    k6 = h * rpf(t+h/2.0,v,r-8.0*k1/27.0+2.0*k2-3544.0*k3/2565.0+
                 1859.0*k4/4104.0-11.0*k5/40.0);
    r1 = r + (25.0*k1/216.0+1408.0*k3/2565.0+2197.0*k4/4104.0-k5/5.0);
    rp =     (16.0*k1/135.0+6656.0*k3/12825.0+28561.0*k4/56430.0-
              9.0*k5/50.0+2.0*k6/55.0);
    rp = rp/h;
    r2 = r + h * rp;

    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 fitz2.c */