restart;
# Determine all two-stage IRK methods of order greater or equal 3
with(LinearAlgebra):
m:=4;

# Since order 3 then can check only for autonomous equation
f:=(t,y)->g(y);

A:=Matrix([[a11,a12],[a21,a22]]);

b:=Vector([b1,b2]);

c:=Vector([a11+a12,a21+a22]);

n:=Dimension(b);

onesvector:=Vector(n,1);
fxi:=Vector(n):
for j from 1 to n do
    fxi[j]:=f(t+c[j]*h,xi[j]);
od:
fxi;
xirhs:=y(t)*onesvector+h*Multiply(A,fxi);



T:=y(t+h)-y(t)-h*Multiply(Transpose(b),fxi);
for j from 1 to m-1
do
    T:=subs(xi=xirhs,T);
od:

T;

S:=series(T,h,m);

eq[1]:=D(t->y(t))(t)=f(t,y(t));

for j from 1 to m-2
do
    eq[j+1]:=simplify(subs(seq(eq[i],i=1..j),D(unapply(eq[j],t))(t)));
od;


T:=simplify(subs(seq(eq[i],i=1..m-1),S));

T1:=coeff(T,h,1);
eq1:=simplify(T1/g(y(t)))=0;


T2:=coeff(T,h,2);
eq2:=simplify(T2/D(g)(y(t))/g(y(t)))=0;


T3:=coeff(T,h,3);

T3a:=coeff(T3,g(y(t)),1);
eq3:=simplify(T3a/D(g)(y(t))^2)=0;


T3b:=coeff(T3,g(y(t)),2);
eq4:=simplify(2*T3b/D(D(g))(y(t)))=0;



solve({eq1,eq2,eq3,eq4},{a11,a12,a21,a22,b1,b2});