(* Code for Exercise 14, Problem 4, multi-variate HE2(1) *)

{a, b} = {0, 100}; (* time interval *)
f[x_, {y_, v_}] := {v, 0.2 (1 - y^2) v - y}; (* r.h.s. function van der Pol *)

Clear[res]; nn = 1000;

Module[{x, y, err, h, hneu, k1, k2, s1, s2, n, print, ag, pg},
 
 (* Initial data *)
 x = 0.; y = {1, -1};
 s1 = 1.; s2 = 1.;
 h = 0.001;
 n = 1;
 res = {};
 ag = 1; pg = 2; (* accuracy and precision goal *)
 
 (* No formatted heads for tabulated data *)
 
 (* Scheme *)
 While[x <= b && n < nn,
  k1 = f[x, y];
  k2 = f[x + h, y + h k1];
  
  (* Heun-Euler 2(1) error formula *)
  err = h ((1/2 - 1) k1 + (1/2 - 0) k2);
  
  hneu = h s1 (Norm[{err[[1]]/(s2 (10^-ag + 10^-pg Abs[y[[1]]])), 
       err[[2]]/(s2 (10^-ag + 10^-pg Abs[y[[2]]]))}, Infinity])^(-1/2);
  
  (* Proceed if err > s2*tol, else reject and repeat with decreased step-size *)
  
  If[ (Norm[{err[[1]]/(s2 (10^-ag + 10^-pg Abs[y[[1]]])), err[[2]]/(
       s2 (10^-ag + 10^-pg Abs[y[[2]]]))}, Infinity])^(-1/2) >= 1.,
   print = "Proceed";
   res = 
    Append[res, {x, y (* StringForm["(``)",Column[y,", "]] *), h, 
      StringForm["(``)", 
       Column[err, 
        ", "]], (Norm[{err[[1]]/(s2 (10^-ag + 10^-pg Abs[y[[1]]])), 
         err[[2]]/(s2 (10^-ag + 10^-pg Abs[y[[2]]]))}, Infinity])^(-1/
       2), hneu, print, StringForm["(``)", Column[k1, ", "]], 
      StringForm["(``)", Column[k2, ", "]]}] ; 
     
   x = x + h;
   y = y + h (1/2) k1 + h (1/2) k2,
   print = "Reject";
   res = 
    Append[res, {x, y (* StringForm["(``)",Column[y,", "]] *), h, 
      StringForm["(``)", 
       Column[err, 
        ", "]], (Norm[{err[[1]]/(s2 (10^-ag + 10^-pg Abs[y[[1]]])), 
         err[[2]]/(s2 (10^-ag + 10^-pg Abs[y[[2]]]))}, Infinity])^(-1/
       2), hneu, print, StringForm["(``)", Column[k1, ", "]], 
      StringForm["(``)", Column[k2, ", "]]}] ; 
  
   ]; 
  
  h = hneu;
  n++;
  ]
 ]