% REFERENCES:
% [1] "Symplectic Integrators" from the MATLAB File Exchange
%      http://www.mathworks.com/matlabcentral/fileexchange/7686

clear;
clf;
cla;
hold off;


n=1;
t_max = 10000;
time_steps=1000000;


% GENERAL SETUP
q0=1.6;
p0=0.5;

% Time interval
tspan = linspace(0,t_max,time_steps);

% SYMPLECTIC SOLVER

% Setup

dpdt=@(q)(-sin(q));

dqdt=@(p)(p);

fprintf('Symplectic Euler (seep):\n');
tic,
[se.q,se.p]=seep(dqdt,dpdt,tspan,q0,p0); % seep is in the package [1]
toc


% EXPLICIT (NON-SYMPLECTIC) Euler
[ee.q,ee.p]=ee(dqdt,dpdt,tspan,q0,p0);

% INBUILT SOLVER
fprintf('ode45:\n');

% Setup
odefun=@(t,y)([y(2);-sin(y(1))]);
y0=[q0;p0];

% Solver options
opt=odeset('Stats','on');

% Solve
tic,
[~,y]=ode45(odefun,tspan,y0,opt); 
toc 


% Convert to p,q format
% evens=2*(1:n);
% odds=2*(1:n)-1;
% ns.q=y(:,odds);
% ns.p=y(:,evens);
ode45ns.q=y(:,1);
ode45ns.p=y(:,2);

 
% PLOT
set(gca,'nextplot','replacechildren'); % axis settings
lsh=plot(ode45ns.q,ode45ns.p,se.q,se.p,'Color',[0 0 1]); % plot
% axis settings
set(lsh(1),'Color',[1 0 0]);
set(lsh(2),'Color',[0 0 1]);
axis square;
axis([-2 2 -2 2]);
title('Phase space evolution; y''''+sin(y) = 0')
xlabel(gca,'q (= y)');
ylabel('p (= y'')');

text('Position',[0.5,0], 'String','Symplectic Euler','Color',[0 0 1]);
text('Position',[0,0], 'String','ode45','Color',[1 0 0]);
text(1.5,-1.9,sprintf('t:%0.1f',tspan(end)));


