% solves Duffing oscillator with Euler

% will need to pass difference between approximations 
% to other subroutines
global diff; 

% allocate memory to the arrays to speed up the code
ntmax=100000;
% sol stores position of particle for all times
sol=zeros(ntmax,1);
energy=zeros(ntmax,1);
% for variable time step need to keep track of deltat and t separately
deltat=zeros(ntmax,1);
t=zeros(ntmax,1);

% measure cputime used: initialize 
cputime1=cputime

% parameters used
tmax=100;
time=0;
dt=0.005;
nt=1;

% initial condition
y=[1.414214;0];
sol(1)=y(1);
t(1)=time;
deltat(1)=dt;

% initial value of total energy
energy1=-y(1)^2/2+y(1)^4/4+y(2)^2/2
% energy(nt) is the deviation of the energy from the initial value
energy(1)=0;

% main computational loop
while time < tmax
 
timep=time+dt;

% to make sure the final time step hits exactly tmax
if timep > tmax
  dt=tmax-time;
end

  [time,y,dt]=stepeuler(time,y,dt);
   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;
 
end
 
time
denergy=energy(nt)
ufinal=sol(nt)

% compute cpu-time used
cputimetotal=cputime-cputime1

% plot the position
figure(1);
plot(t(1:nt),sol(1:nt));
xlabel('time');
ylabel('u');
% plot of time steps 
figure(2);
plot(t(1:nt),deltat(1:nt),'+');
xlabel('time');
ylabel('delta t');
% plot of change in energy (for monitoring quality)
figure(3);
plot(t(1:nt),energy(1:nt));
xlabel('time');
ylabel('Energy');
 

