func real Y_(T, Y);
real T, Y;
return T*sqrt(Y);

def  DT = 0.1;
real T, Y, Exact, DY1, DY2, DY3, DY4;
[Text(0, "  T       RK           Exact        Error^m^j");
T:= 0.;  Y:= 1.;
repeat  if Mod(T+.001, 1.) < .01 then
            [Format(2, 1);
            RlOut(0, T);
            Format(5, 7);
            RlOut(0, Y);
            Exact:= sq(T*T+4.)/16.;
            RlOut(0, Exact);
            RlOut(0, Y-Exact);
            CrLf(0);
            ];
        DY1:= DT * Y_(T, Y);
        DY2:= DT * Y_(T+DT/2., Y+DY1/2.);
        DY3:= DT * Y_(T+DT/2., Y+DY2/2.);
        DY4:= DT * Y_(T+DT, Y+DY3);
        Y:= Y + (DY1 + 2.*DY2 + 2.*DY3 + DY4) / 6.;
        T:= T + DT;
until   T > 10.;
]