Runge_Kutta: procedure options (main);           /* 10 March 2014 */
   declare (y, dy1, dy2, dy3, dy4) float (18);
   declare t fixed decimal (10,1);
   declare dt float (18) static initial (0.1);

   y = 1;
   do t = 0 to 10 by 0.1;
      dy1 = dt * ydash(t, y);
      dy2 = dt * ydash(t + dt/2, y + dy1/2);
      dy3 = dt * ydash(t + dt/2, y + dy2/2);
      dy4 = dt * ydash(t + dt,   y + dy3);

      if mod(t, 1.0) = 0 then
         put skip edit('y(', trim(t), ')=', y, ', error = ', abs(y - (t**2 + 4)**2 / 16 ))
                      (3 a, column(9), f(16,10), a, f(13,10));
      y = y + (dy1 + 2*dy2 + 2*dy3 + dy4)/6;
   end;


ydash: procedure (t, y) returns (float(18));
   declare (t, y) float (18) nonassignable;
   return ( t*sqrt(y) );
end ydash;

end Runge_kutta;