% Commutator free method of order 4

tend=5;
epsilon=0.05;
omega=1/epsilon;
h=0.01;
y0=[1;0];
Nsteps = round(tend/h);
Y=zeros(Nsteps,2);
Y(1,:)=y0';
oA = omega*[0,1;-1,0];

for j=1:Nsteps
  y1 = cfree4step(Y(j,:)',h,oA,epsilon);
  Y(j+1,:)=y1';
end
figure(1)
plot(Y(:,1),Y(:,2),'k-');
title('Phase plot with CFREE4 method');
    
tend=1;

func=@(t,y) [omega*y(2)+4*y(2)^3-epsilon*y(1);-omega*y(1)-4*y(1)^3-epsilon*y(2)];
options = odeset('RelTol',1e-10,'AbsTol',1e-11);
[Tex,Yex] = ode45(func,[0,tend],y0,options);
HH=zeros(6,1); ERR=zeros(6,1);
for k=0:5,
    Nsteps=10*2^k;
    h=tend/Nsteps;
    Yc=y0;
    for j=1:Nsteps,
      Yc=cfree4step(Yc,h,oA,epsilon);
    end
    HH(k+1,1)=h;
    ERR(k+1,1)=norm(Yc-Yex(end,:)',2);
end
figure(2)
loglog(HH,ERR,'b-o',HH,HH.^4,'r--');
legend('Glob err at t=1','Ref line order 4','Location','North')