restart;
k1 := f(x, y(x)); k2 := f(x+a1*h, y(x)+b1*h*k1); k3 := f(x+a2*h, y(x)+b2*h*k2);
`output redirected...`print(); # input placeholder
T := y(x+h)-y(x)-h*(g1*k1+g2*k2+g3*k3);
`output redirected...`print(); # input placeholder
R1 := series(T, h, 5);
s1 := diff(y(x), x) = (D(y))(x);
t1 := diff(y(x), x) = f(x, y(x));
`output redirected...`print(); # input placeholder
eq[1] := subs(s1, t1);
`output redirected...`print(); # input placeholder
for i from 2 to 5 do t1 := diff(eq[i-1], x); eq[i] := subs(s1, t1) end do;
??print(); # input placeholder
eq[1];
`output redirected...`print(); # input placeholder
eq[2];
`output redirected...`print(); # input placeholder
eq[3];
`output redirected...`print(); # input placeholder
s2 := [seq(eq[6-j], j = 1 .. 5)];
??print(); # input placeholder
# 
# 
test2 := subs(s2, coeff(P1, h, 3));
P1 := subs(s2, subs(s2, subs(s2, subs(s2, convert(R1, polynom)))));
??print(); # input placeholder
coeff(P1, h, 0);
`output redirected...`print(); # input placeholder
t2 := coeff(P1, h, 1);
`output redirected...`print(); # input placeholder
e1 := coeff(t2, f(x, y(x)), 1);
`output redirected...`print(); # input placeholder
S1 := solve(e1, {g1});
`output redirected...`print(); # input placeholder
P2 := simplify(subs(S1, P1));
??print(); # input placeholder
coeff(P2, h, 0);
`output redirected...`print(); # input placeholder
coeff(P2, h, 1);
`output redirected...`print(); # input placeholder
t3 := coeff(P2, h, 2);
`output redirected...`print(); # input placeholder
e2 := coeff(t3, (D[1](f))(x, y(x)));
`output redirected...`print(); # input placeholder
S2 := solve(e2, {a1});
`output redirected...`print(); # input placeholder
P3 := simplify(subs(S2, P2));
??print(); # input placeholder
t4 := coeff(P3, h, 2);
`output redirected...`print(); # input placeholder
e3 := simplify(coeff(t4, (D[2](f))(x, y(x)))/f(x, y(x)));
`output redirected...`print(); # input placeholder
S3 := solve(e3 = 0, {b1});
`output redirected...`print(); # input placeholder
P4 := simplify(subs(S3, P3));
??print(); # input placeholder
coeff(P4, h, 2);
`output redirected...`print(); # input placeholder
S1, S2, S3;
`output redirected...`print(); # input placeholder
t5 := coeff(P4, h, 3);
`output redirected...`print(); # input placeholder
e4 := simplify(coeff(t5, (D[1](f))(x, y(x)))/(D[2](f))(x, y(x)));
`output redirected...`print(); # input placeholder
S4 := solve(e4 = 0, {a2});
`output redirected...`print(); # input placeholder
P5 := simplify(subs(S4, P4));
t6 := coeff(P5, h, 3);
`output redirected...`print(); # input placeholder
e5 := simplify(coeff(t6, (D[2](f))(x, y(x)), 2)/f(x, y(x)));
`output redirected...`print(); # input placeholder
S5 := solve(e5 = 0, {g2});
`output redirected...`print(); # input placeholder
P6 := simplify(subs(S5, P5));
??print(); # input placeholder
t7 := coeff(P6, h, 3);
`output redirected...`print(); # input placeholder
e6 := coeff(t7, (D[1, 1](f))(x, y(x)));
`output redirected...`print(); # input placeholder
S6 := solve(e6 = 0, {g3});
`output redirected...`print(); # input placeholder
S6b := S6[2];
`output redirected...`print(); # input placeholder
P7 := simplify(subs(S6b, P6));
??print(); # input placeholder
coeff(P7, h, 3);
`output redirected...`print(); # input placeholder
S1, S2, S3, S4, S5, S6;
`output redirected...`print(); # input placeholder
coeff(P7, h, 4);
??print(); # input placeholder
#