%%
%----------------------------%
% Example Runge-Kutta method %
%----------------------------%
%----------------%
% Exercise 6.0.2 %
%----------------%
clear all
close all
clc

a=1; b=4; M=5;
t=linspace(a,b,M);
dt=t(2)-t(1);

f=@(t)(-2*(t-5/2));
x1=-9/4;


xt_E(1)=x1;
xt_modE(1)=x1;
xt_RK(1)=x1;

for i=1:M-1
   xt_E(i+1)   =xt_E(i) + dt*f(t(i));
   xt_modE(i+1)=xt_modE(i) + dt*( f(t(i)) + f(t(i+1)) )/2;
   xt_RK(i+1)  =xt_RK(i) + dt*f(t(i)+dt/2);
end

N=100;
tN=linspace(a,b,N);
xt_exact=@(t)( -(t-5/2).^2 );

hold on
plot(tN,xt_exact(tN),'r-','LineWidth',2)
plot(t,xt_modE,'c-','Linewidth',3)
plot(t,xt_RK,'k--')
plot(t,xt_E,'b-')
legend('exact solution','2nd order Runge-Kutta', 'modified Euler', '1st order Euler',...
    'Location','South')
xlabel('t','fontsize',16),ylabel('x(t)','fontsize',16)
hold off

print -f -depsc 'RK_approx.eps'
%%
xt_ex=-(t-5/2).^2;

max(abs(xt_ex-xt_E))
dt
max(abs(xt_ex-xt_RK))
(dt)^2
%%
%-----------------------------%
% Exercise 6.0.3              %
% now with external functions %
%-----------------------------%
clear all
close all
clc

a=1; b=4; M=5;
t=linspace(a,b,M);
dt=t(2)-t(1);

f=@(x,t)(-2*(t-5/2));
x1=-9/4;

[xt_RK, t] = RK_method(a, x1, b, M, f);

N=100;
tN=linspace(a,b,N);
xt_exact=@(t)( -(t-5/2).^2 );

hold on
plot(tN,xt_exact(tN),'r-','LineWidth',2)
%plot(t,xt_modE,'c-','Linewidth',3)
plot(t,xt_RK,'k--')
%plot(t,xt_E,'b-')
legend('exact solution','2nd order Runge-Kutta', 'modified Euler', '1st order Euler',...
    'Location','South')
xlabel('t','fontsize',16),ylabel('x(t)','fontsize',16)
hold off

%print -f -depsc 'RK_approx.eps'
%%
clear all
close all
clc

t1=0; tN=1; N=7;
x1=1;

t=linspace(t1,tN,N);
dt=t(2)-t(1);

f=@(x)(-3*x);

xt_E(1)=x1;
xt_RK(1)=x1;
xt_modE(1)=x1;
for i=1:N-1
   xt_E(i+1)   =xt_E(i) + dt*f(xt_E(i));
   
   xt_intE(i)=xt_modE(i) + dt*f(xt_modE(i))
   xt_modE(i+1)=xt_modE(i) + dt*( f(xt_modE(i))+f(xt_intE(i)) )/2;
   
   xmid(i)     =xt_RK(i)+(dt/2)*f(xt_RK(i));
   xt_RK(i+1)  =xt_RK(i) + dt*f(xmid(i));
end

Npl=100; tpl=linspace(t1,tN,Npl);

hold on
plot(tpl,exp(-3*tpl),'r-','LineWidth',2)
plot(t,xt_E,'b-')
plot(t,xt_modE,'c-','Linewidth',3)
plot(t,xt_RK,'k--')
legend('exact solution', 'Euler approximation',...
    'modified Euler approximation', '2nd order Runge-Kutta', 'Location','North')
xlabel('t','fontsize',16),ylabel('x(t)','fontsize',16)
hold off

print -f -depsc 'exp-decay_RK_approx.eps'
%%
%-----------------------------%
% Exercise 6.0.4              %
% now with external functions %
%-----------------------------%
clear all
close all
clc

t1=0; tN=5; N=20;
x1=1;
f=@(x,t)(5*x-x.^2);

[xt_RK, t] = RK_method(t1, x1, tN, N, f);
[xt_E, t] = Euler_method(t1, x1, tN, N, f);
[xt_modE, t] = modEuler_method(t1, x1, tN, N, f);

N=100;
tN=linspace(t1,tN,N);
xt_exact=@(t)( 5*exp(5*t)./(4+exp(5*t)) );

hold on
plot(tN,xt_exact(tN),'r-','LineWidth',2)
plot(t,xt_RK,'k--')
plot(t,xt_modE,'c-','Linewidth',3)
plot(t,xt_E,'b-')
legend('exact solution','2nd order Runge-Kutta', 'modified Euler', '1st order Euler',...
    'Location','South')
xlabel('t','fontsize',16),ylabel('x(t)','fontsize',16)
hold off


%%
%-----------------------------%
% Exercise 7.0.1              %
%-----------------------------%
clear all
close all
clc

t1=0; tN=4*pi; N=50;
t=linspace(t1,tN,N); dt=t(2)-t(1);
x0=1; v0=0; x_exact=@(t)( cos(t) );

x(1)=x0; v(1)=v0;
for i=1:N-1
    xm(i)=x(i) + (dt/2)*v(i);
    vm(i)=v(i) + (dt/2)*(-x(i));
    
    x(i+1)=x(i) + dt*vm(i);
    v(i+1)=v(i) + dt*(-xm(i)); 
end

figure(1)
hold on
plot(t,x_exact(t),'r-','LineWidth',2)
plot(t,x,'g-')
legend('exact solution',...
    'Runge-Kutta approximation',...
    'Location','NorthWest')
hold off





