% fitz.m  FitzHugh-Nagumo equations  system of two ODE's
function fitz
  h = 0.001;
  t = 0.0;
  step = 0.03125;
  td = step;
  v = -1.0;
  r = 1.0;
  j = 1;
  while t<25
    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(abs(v1-v2) + abs(r1-r2) > 0.000001)
      h = h/2.0;
      sprintf('h reduced to %g \n', h)
      if(h<1.0e-6)
        break
      end
      continue
    end
    v = v2;
    r = r2;
    t = t + h;
    if(t>=td)
      tt(j) = t;
      vv(j) = v;
      rr(j) = r;
      j = j + 1;
      td = td + step;
    end
  end
  axis tight
  figure(1)
  plot(tt, vv, '-')
  axis([0 25 -12.5 12.5])
  title('v vs t')
  xlabel('t')
  ylabel('v')
  figure(2)
  plot(tt, rr, '-')
  axis([0 25 -12.5 12.5])
  title('r vs t')
  xlabel('t')
  ylabel('r')
  return
  function z = vpf(v, r)
    z = v - v*v*v/3.0 + r;
  return
  end
  function z = rpf(v, r)
    z = -(v - 0.2 - 0.2*r);
  return
  end
end