% solves Duffing oscillator with Euler or RK4 with
% time-step control and Richardson extrapolation

format long
flops1=flops

global diff;

tmax=100;
y=[1.414214;0];
time=0;
tol=0.15e-7;
dt=0.05;
nt=1;
sol(1)=y(1);
t(1)=time;
deltat(1)=dt;
energy1=-y(1)^2/2+y(1)^4/4
energy(1)=0;
difference(1)=0;

while time < tmax
 
timep=time+dt;

if timep > tmax
  dt
  time
  dt=tmax-time
end

   [time,y,dt,success]=steprk(time,y,dt,tol);
%  [time,y,dt,success]=steprkoc(time,y,dt,tol);
%   [time,y,dt,success]=steprkocwe(time,y,dt,tol);
%  [time,y,dt,success]=stepeuleroc(time,y,dt,tol);
  if success > 0
   nt=nt+1;
   sol(nt)=y(1);
   t(nt)=time;
   energy(nt)=-y(1)^2/2+y(1)^4/4+y(2)^2/2-energy1;
   deltat(nt)=dt;
   difference(nt)=diff;
  end

end
 
nt
time
 
denergy=energy(nt)
ufinal=sol(nt)
deltat2=deltat(floor(nt/4):floor(3*nt/4));
dtmax=max(deltat2)
dtmin=min(deltat2)

flops2=flops;
flop=flops2-flops1
flops1=flops2;

figure(1);
plot(t(1:nt),sol(1:nt));
ylabel('u');
figure(2);
plot(t(1:nt),deltat(1:nt),'+');
ylabel('Delta t');
figure(3);
plot(t(1:nt),energy(1:nt));
ylabel('Energy');
figure(4);
plot(t(1:nt),difference(1:nt));
ylabel('Difference in u');
