for i from 1 to 2 do
  k[i]:=f(y(t0)+h*sum(a[i,j]*K[j],j=1..2));
od;
for i from 1 to 2 do
  kn[i]:=subs({K[1]=k[1],K[2]=k[2]},k[i]);
od;
for i from 1 to 2 do
  k[i]:=subs({K[1]=kn[1],K[2]=kn[2]},kn[i]);
od;
D(y):=t->f(y(t));
y1:=y(t0)+h*sum(b[j]*k[j],j=1..2);;
p:=collect(convert(taylor(y(t0+h)-y1,h,5),polynom),h);
vars:={a[1,1],a[1,2],a[2,1],a[2,2],b[1],b[2]};
eqns:={coeffs(expand(p),indets(p) minus vars)};
sols:=solve(eqns,vars);
allvalues(sols);