(* Code for HE2(1) *)

{a, b} = {0., 8}; (* interval on x-axis *)
f[x_, y_] := 1. (1. - 1./4 Cos[y])^2; (* r.h.s. function of the differential equation *)

Clear[res]; nn = 200;

Module[{x, y, err, h, hneu, k1, k2, s1, s2, n, print, ag, pg},
 
 (* Initial data *)
 x = 0.; y = 0.;
 s1 = 1.; s2 = 1.;
 h = 0.001;
 n = 1;
 ag = 4; pg = 4;
 
 (* no row with formatted titles *)
  
 (* Scheme *)
 While[x <= b + h && 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 (Abs[(err/(s2 (10^-ag + 10^-pg Abs[y])))])^(-1/2);
  
  (* Proceed if err > s2*tol, else reject and repeat with decreased step-size *)
  
  If[ (Abs[(err/(s2 (10^-ag + 10^-pg Abs[y])))])^(-1/2) >= 1.,
   print = "Proceed";
   res = Append[
     res, {x, y, h, 
      err, (Abs[(err/(s2 (10^-ag + 10^-pg Abs[y])))])^(-1/2), hneu, 
      print, k1, k2}] ;
   x = x + h;
   y = y + h (1/2) k1 + h (1/2) k2,
   
   print = "Reject";
   res = Append[
     res, {x, y, h, 
      err, (Abs[(err/(s2 (10^-ag + 10^-pg Abs[y])))])^(-1/2), hneu, 
      print, k1, k2}] ;
   ];
  
  h = hneu;
  n++;
  ]
 ]