Question: solution of the differential equation by the Milne method

Good evening!!! I have a task to implement the task of Cauchy by the method of Milne, wrote the code, but did not understand it until the end, help to understand? what's wrong?
First calculate four "initial" values by the method of Runge-Kutta methods, then use the method of Milne, the Fact that two times running, perhaps extra?

restart;
with(plots):
a:=0; b:=1; eps:=evalf(10^(-3)):
f:=unapply(2*x*(x^2+y),x,y);
G:=simplify(dsolve({diff(y(x),x)=f(x,y(x)),y(a)=1}));                    
N:=15: h:=(b-a)/N:
for i from 0 to N do 
x[i]:=a+i*h: 
end do:
y[0]:=1;
s[0]:=1;
for i from 0 to 2 do 
t[1]:=evalf(h*f(x[i],y[i])):
t[2]:=evalf(h*f(x[i]+h/2,y[i]+t[1]/2)): 
t[3]:=evalf(h*f(x[i]+h/2,y[i]+t[2]/2)):
t[4]:=evalf(h*f(x[i]+h,y[i]+t[3])):
y[i+1]:=evalf(y[i]+(t[1]+2*t[2]+2*t[3]+t[4])/6):
q[1]:=evalf(h*f(x[i],s[i])):
q[2]:=evalf(h*f(x[i]+h/2,s[i]+q[1]/2)): 
q[3]:=evalf(h*f(x[i]+h/2,s[i]+q[2]/2)):
q[4]:=evalf(h*f(x[i]+h,s[i]+q[3])):
s[i+1]:=evalf(s[i]+(q[1]+2*q[2]+2*q[3]+q[4])/6):
end do;
for i from 3 to N-1 do 
y[i+1]:=evalf(y[i-3]+((4*h)/3)*(2*f(x[i],y[i])-f(x[i-1],y[i-1])+2*f(x[i-2],y[i-2]))):
s[i+1]:=evalf(s[i-1]+(h/3)*(f(x[i+1],y[i+1])+4*f(x[i],s[i])+f(x[i-1],s[i-1]))):
d[i+1]:=abs(y[i+1]-s[i+1])/29:
if abs(d[i+1]) < eps then y[i]:=y[i]:
else y[i]:=s[i];
end if: end do;
s1:=plot(rhs(G),x=a..b,color=yellow):
s2:=pointplot({seq([x[k],y[k]],k=0..N)}): 
display(s1,s2);