function [tb,u,v] = RK4System_oscillator(int,u0,v0,nsteps)
% Compares the Forward Euler, Runge-Kutta and Symplectic Euler with the exact solution
% (for a system)
% when
% solving the simple oscillator problem : u"+u = 0, u(0)=-1,u'(0)=-2,
% as a system of first-order equations: u'=v, v'=-u, u(0)=-1, v(0)=-2.
% The exact solution is: u(t)= -2*sin(t)-cos(t); v(t) = -2cos(t)+sin(t).
%
% Input: Interval [ta,T], initial condition x0,  and number of steps n
% Output: the approximate solution u,v, at time T
%
% Example usage: [T,u,v] = RK4System_oscillator([0 2*pi],-1,-2,floor(2*pi/0.1))
%
% Suggestion: to observe convergence, run with smaller and smaller
% time-steps, e.g.: RK4System_oscillator([0 2*pi],-1,-2,2*floor(2*pi/0.1))
% and RK4System_oscillator([0 2*pi],-1,-2,4*floor(2*pi/0.1))
%
% Remark: run on the trible the interval to see how 
% the numerical solutions depart from the exact solution, e.g., 
% RK4System_oscillator([0 6*pi],-1,-2,floor(6*pi/0.1)+1)
% !!!!!!!!!!!!!!!!!!!!!!!!
% Notice that the first order Symplectic Euler is more accurate than the
% fourth order RK4. Why??? 
% !!!!!!!!!!!!!!!!!!!!!!!!


ta = int(1); tb = int(2);
h = (tb - ta)/nsteps; %%step size h=(T-ta)/n
time=zeros(nsteps,1);
%Initial conditions

%RK4
u(1) = u0;
v(1) = v0;
%forward Euler
x(1) = u0;
y(1) = v0;
%symplectic Euler
seu(1) = u0;
sev(1) = v0;
%exact solution
uexact(1) = u0;
vexact(1) = v0;

time(1) = ta;


format long 
%RK4 solution
for j = 1:nsteps
    K11 = h*xdot(time(j),u(j),v(j));
    K12 = h*ydot(time(j),u(j),v(j));
    
    K21 = h*xdot(time(j)+h/2,u(j)+K11/2,v(j)+K12/2);
    K22 = h*ydot(time(j)+h/2,x(j)+K11/2,y(j)+K12/2);
    
    K31 = h*xdot(time(j)+h/2,u(j)+K21/2,v(j)+K22/2);
    K32 = h*ydot(time(j)+h/2,u(j)+K21/2,v(j)+K22/2);
    
    K41 = h*xdot(time(j)+h,u(j)+K31,v(j)+K32);
    K42 = h*ydot(time(j)+h,u(j)+K31,v(j)+K32);
    
    u(j+1) = u(j) + (K11 + 2*K21 + 2*K31 + K41)./6;
    v(j+1) = v(j) + (K12 + 2*K22 + 2*K32 + K42)./6;

%forward Euler solution
    x(j+1) = x(j) + h*xdot(time(j),x(j),y(j));
    y(j+1) = y(j) + h*ydot(time(j),x(j),y(j));
%symplectic Euler solution
    seu(j+1) = seu(j) + h*xdot(time(j),seu(j),sev(j));
    sev(j+1) = sev(j) + h*ydot(time(j),seu(j+1),sev(j));
%exact solution
    uexact(j) = -2*sin(time(j))-cos(time(j));
    vexact(j) = -2*cos(time(j))+sin(time(j));

    time(j+1) = time(j) + h;
end
    uexact(nsteps+1) = -2*sin(time(nsteps+1))-cos(time(nsteps+1));
    vexact(nsteps+1) = -2*cos(time(nsteps+1))+sin(time(nsteps+1));

figure(1)
    title('System of 2 equations with RK4')
    hold all
    plot(time,u,'.-.','LineWidth',1,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor','g',...
                'MarkerSize',10);

    plot(time,v,'-','LineWidth',2,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor','g',...
                'MarkerSize',10);
    plot(time,x,'.-.','LineWidth',1,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor','g',...
                'MarkerSize',10);

    plot(time,y,'-','LineWidth',2,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor','g',...
                'MarkerSize',10);
    plot(time,seu,'.-.','LineWidth',1,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor','g',...
                'MarkerSize',10);

    plot(time,sev,'-','LineWidth',2,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor','g',...
                'MarkerSize',10);        
    plot(time,uexact);

    plot(time,vexact);
legend('RK4-u','RK4-v','Euler-x', 'Euler-y','Symplectic Euler u','Symplectic Euler v','uexact','vexact')

figure(2)
hold all
plot(u,v,'.-.') % RK$ solution
plot(x,y,'o')
plot(seu,sev)
plot(uexact,vexact)
legend('RK4:u-v', 'Euler:x-y','symplectic Euler','exact')
axis on
axis tight
title('Phase plane u-v')


for i=1:nsteps
    figure(3)
    hold all
%    plot(uexact(i),vexact(i),'.-.') %exact
    plot(u(i),v(i),'.') %RK4
    plot(x(i),y(i),'o')  %Forward Euler
    plot(seu(i),sev(i),'*') % symplectic Euler

    axis([-4.5 4.5 -4.5 4.5])
    pause
    
%    legend('exact','RK4:u-v', 'Euler:x-y','symplectic Euler')
end
end



function z1 = xdot(t,x,y) %the RHS of the ODE
z1 = y;
end

function z2 = ydot(t,x,y)
z2 = -x;
end
