restart;
# Improved more Efficient Version
with(LinearAlgebra):
m:=5;

A:=Matrix([[0,0],[1/3,1/3]]);

b:=Vector([1/4,3/4]);

c:=Vector([0,2/3]);

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;
xirhs:=y(t)*onesvector+h*Multiply(A,fxi);





xilhs:=Vector(n,i->xi[i]):
for j from 1 to m-1
do
    xilhs:=subs(xi=xilhs,xirhs);
    xilhs:=Map(x->series(x,h,j),xilhs);
    xilhs:=Map(x->convert(x,polynom),xilhs);
od:
xilhs;

R:=Vector(subs(xi=xilhs,fxi));

T:=y(t+h)-y(t)-h*Multiply(Transpose(b),R);

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;



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