function y1 = cfree4step(y0,h,oA,epsilon)
%
%  Take one step with the cfree4 method (by Celledoni, Marthinsen and
%  Owren)
%
Y1=y0;
k1=G(Y1,epsilon);
Y2=affexp(oA,k1,Y1,h/2);
k2 = G(Y2,epsilon);
Y3 = affexp(oA,k2,Y1,h/2);
k3 = G(Y3,epsilon);
Y4 = affexp(oA/2,k3-k1/2,Y2,h);
k4 = G(Y4,epsilon);
yp = affexp(oA/2,(3*k1+2*k2+2*k3-k4)/12,Y1,h);
y1 = affexp(oA/2,(-k1+2*k2+2*k3+3*k4)/12,yp,h);

end

function Fr = G(y,epsilon)
Fr=[4*y(2)^3-y(1)*epsilon; -4*y(1)^3-y(2)*epsilon];
end

function y1=affexp(A,b,y0,h)
% In this case we know that A is nonsingular, but in general, take care if
% A has tiny eigenvalues
eA=expm(h*A);
y1 = eA*y0 + (eA-eye(2))*(A\b);
end