restart:

N := 20: RKNOrdnung:= proc(f1, f2, x0, y0, a, b, n) local j,h,t,x,y,k1x,k2x,k3x,k4x,k1y,k2y,k3y,k4y; h:=(b-a)/n; t[0]:=a: x[0]:=x0: y[0]:=y0: for j from 0 by 1 while j<n do k1x:=evalf( h*f1(t[j], x[j], y[j] ) ); k1y:=evalf( h*f2(t[j], x[j], y[j] ) ); k2x:=evalf( h*f1(t[j]+h/2, x[j]+k1x/2, y[j]+k1y/2 ) ); k2y:=evalf( h*f2(t[j]+h/2, x[j]+k1x/2, y[j]+k1y/2 ) ); k3x:=evalf( h*f1(t[j]+h/2, x[j]+k2x/2, y[j]+k2y/2 ) ); k3y:=evalf( h*f2(t[j]+h/2, x[j]+k2x/2, y[j]+k2y/2 ) ); k4x:=evalf( h*f1(t[j], x[j]+k3x, y[j]+k3y) ); k4y:=evalf( h*f2(t[j], x[j]+k3x, y[j]+k3y) ); t[j+1]:=t[j]+h; x[j+1]:=x[j]+1.0/6.0*(k1x+2*k2x+2*k3x+k4x); y[j+1]:=y[j]+1.0/6.0*(k1y+2*k2y+2*k3y+k4y); end do: t,x,y; end proc:

f1:=(x,y,yStrich) -> yStrich; f2:=(x,y,yStrich) -> 1/8*(32+2*x^3-y*yStrich);

NiM+SSNmMUc2ImYqNiVJInhHRiVJInlHRiVJKHlTdHJpY2hHRiVGJTYkSSlvcGVyYXRvckdGJUkmYXJyb3dHRiVGJTkmRiVGJUYl

NiM+SSNmMkc2ImYqNiVJInhHRiVJInlHRiVJKHlTdHJpY2hHRiVGJTYkSSlvcGVyYXRvckdGJUkmYXJyb3dHRiVGJSwoIiIlIiIiKiQ5JCIiJCNGMEYvKiY5JUYwOSZGMCMhIiIiIilGJUYlRiU=

alpha:=17: beta:=43/3: t[0]:=alpha: t[1]:=beta:

for i from 2 to 10 do: t1,x1,y1 := RK(f1, f2, alpha, t[i-1], alpha, beta, N); t2,x2,y2 := RK(f1, f2, alpha, t[i-2], alpha, beta, N); t[i]:=t[i-1] - (y1[N]-beta)/(y1[N]-y2[N]); end do: