(* PROJECT: INITIAL VALUE PROBLEMS
EULER' s - IMPROVED EULER' S - RUNGE KUTTA 
MILAN RAMAIYA - 2001.12  .06 *)
  
Eulers[f_,  x0_,  y0_,  xn_,  
            h_]  := 
(
yi = 
              y0; 
For[xi = x0,  
              xi ≤ (xn - h),  
              xi += h,  yi += h*f[xi, yi];]; 
Return[
              yi];
);
  
EulersImp[f_,  x0_,  y0_,  xn_,  
            h_]  := 
(
yi = 
              y0; 
For[xi = x0,  
              xi ≤ (xn - h),  
              xi += h, 
yb = 
                yi + h*f[xi, 
                      yi]; 
yi += (h/2) (f[xi, yi] + 
                      f[xi + h, yb]);]; 
Return[
              yi];
);
  
RKutta[f_,  x0_,  y0_,  xn_,  
          h_]  := 
(
yi = 
            y0; 
For[xi = x0,  
            xi ≤ (xn - h),  
            xi += h,  
K1 = 
              h*f[xi, yi]; 
K2 = 
              h*f[xi + (h/2), yi + (K1/2)]; 
K3 = 
              h*f[xi + (h/2), yi + (K2/2)]; 
K4 = 
              h*f[xi + h, 
                  yi + K3]; 
yi += ((K1 + 2  K2 + 2  K3 + K4)/6);
]; 
Return[
            yi];
);

fp[x_, y_] = y - (x^2) + 1;
Print[
        Eulers[fp,  0,   .5,   .5,   .05]];
Print[
        EulersImp[fp,  
          0,   .5,   .5,   .05]];
Print[RKutta[fp,  0,   .5,   .5,   .05]];