restart;
#
a := -1.;
#
b := 0.;
#
H := 0.786;
t := 1/2;
#
yy:=(x)->sin(x^3+1);
#
y0:=yy(a);
#
F:=(x,y)->3*x^2*cos(x^3+1)-sin(x^3+1)+y;
#
eps:=0.0001;
#
RGKSTEP := proc(XOLD, YOLD, H, t)
local K1, K2, K3, K4, xnew, ynew;
K1 := H * F(XOLD,YOLD);
K2 := H * F(XOLD + H/2, YOLD + K1/2);
K3 := H * F(XOLD + H, YOLD + 2 * K2 - K1);
ynew := evalf(YOLD + (1/6) * (K1 + 4 * K2 + K3));
xnew := XOLD + H;
return (xnew, ynew);
end proc:
#
RGKSTEP(a, y0, H, t);
yy(a+H);
#
RGKH := proc(X0, Y0, B, n, t)
local X, Y, xx, yy, i, H;
H:=(B-X0)/n;
X := array(0..n):
Y := array(0..n):
X[0] := X0;
Y[0] := Y0;
xx := X0;
yy := Y0;
for i from 1 to n do
   X[i], Y[i] := RGKSTEP(xx, yy, H, t);
   xx := X[i];
   yy := Y[i];
end do:
return X, Y;
end proc:
#
N:=5;
#
XX, YY := RGKH(a, y0, b, N, t);
#
for i from 0 to N do
   print("x=", XX[i], " yrk=",YY[i], " yt=",yy(XX[i])," pogr=",abs(YY[i]-yy(XX[i])));
end do:
#
RGKAUTOH := proc(X0, Y0, H, B, EPS, t)
local h, X, Y, xx, yy, y1, pr, y2, xi, k, i, cd, xnew;
X[0] := X0;
Y[0] := Y0;
xx := X0;
yy := Y0;
h := H;
i := 1;
cd:=0;
while (xx + h <= B) and (cd =0) do
  xnew, y1 := RGKSTEP(xx, yy, h, t);
  xnew, pr := RGKSTEP(xx, yy, h/2, t);
  xnew, y2 := RGKSTEP(xx + h/2, pr, h/2, t);
  k := 0;
  while (k < 5) and (abs ((y1 - y2)/7) >= EPS) do
    k := k + 1;
    h := h/2;
    y1:= pr;
    xnew, pr := RGKSTEP(xx, yy, h/2, t);
    xnew, y2 := RGKSTEP(xx + h/2, pr, h/2, t);
  end do:
  if k < 5 then
    X[i] := xx + h;
    Y[i] := y2;
    xx := X[i];
    yy := Y[i];
    i := i + 1;
    if (abs((y1 - y2)/(7)) < EPS/16) then
      h := h*2;
    end if:
    if (xx + h > B) then
      h := B - xx;
      if (h = 0) then
        cd := 99;
      fi:
    fi:
  else
    cd := 99;
  end if:
end do:

return X, Y, i-1 ;
end proc:
#
h0:=0.01;
a; b; yy(a);
#
Xh, Yh, ih := RGKAUTOH(a, yy(a), h0, b, 0.000001, t):
print ("Значение X: ",Xh[0],
       " Посчитанное значение в точке X: ", Yh[0],
       " Точное значение в X: ", evalf(yy(Xh[0])),
       " Шаг: ", evalf(h0),
       " Погрешность: ", abs(Yh[0]- evalf(yy(Xh[0])))
      );
for i from 1 to ih  do
   print ("Значение X: ", Xh[i],
          " Посчитанное значение в точке X: ", Yh[i],
          " Точное значение в X: ", evalf(yy(Xh[i])),
          " Шаг: ", evalf(Xh[i]-Xh[i-1]),
          " Погрешность: ", abs(Yh[i]- evalf(yy(Xh[i])))
         );
end do:
#
with(plots):
gr1:=pointplot({seq([Xh[i], Yh[i]], i=0..ih)}, color = red):
gr2:=plot(yy(x), x=a..b, color=blue, thickness=2):
display(gr1, gr2);