%**************************************************************************
%*     Problem 5.1                                                        *
%*     filename: problem05_1.m                                            *
%*     Programed by Ryoichi Kawai for Computational Physics Course        *
%*     Revised on 02/07/2016.                                             *
%**************************************************************************
clear all;

% parameters

tau=1.0; % time is measured with tau as unit.
g0=0.1;  % initial frictional constant
m=1.0;   % mass

% Control parameters
tmax=10; N=1000; h=tmax/N;

% initial conditions
t(1)=0;
v_rk2(1)=10.0;
v_rk4(1)=10.0;

% array for errors
err_rk2(1)=0.0;
err_rk4(1)=0.0;

% 2nd-order Runge-Kutta method
for n=1:N-1
    t(n+1) = t(1)+n*h;
    k1 = -g0*exp(-t(n))*v_rk2(n);
    v_mid = v_rk2(n)+k1*h/2;
    k2 = -g0*exp(-(t(n)+h/2))*v_mid;
    v_rk2(n+1)=v_rk2(n)+k2*h;
    err_rk2(n+1)=v_rk2(n+1)-v_rk2(1)*exp(-g0/m*(1-exp(-t(n+1))));
end

% 4th-order Runge-Kutta method
for n=1:N-1
    k1 = -g0*exp(-t(n))*v_rk4(n);
    v_mid = v_rk4(n)+k1*h/2;
    k2 = -g0*exp(-(t(n)+h/2))*v_mid;
    v_mid = v_rk4(n)+k2*h/2;
    k3 = -g0*exp(-(t(n)+h/2))*v_mid;
    v_end = v_rk4(n)+k3*h;
    k4 = -g0*exp(-(t(n+1)))*v_end;
    v_rk4(n+1)=v_rk4(n)+ (k1+2*(k2+k3)+k4)*h/6;
    err_rk4(n+1)=v_rk4(n+1)-v_rk4(1)*exp(-g0/m*(1-exp(-t(n+1))));
end

% plot the velocity
subplot(1,2,1); 
p=plot(t,v_rk2,t,v_rk4);
xlabel('t');
ylabel('v');
set(p(1),'Color','Blue','Linewidth',2);
set(p(2),'Color','Red','Linewidth',2);
legend(p,{'RK2','RK4'});

% plot the absolute errors
subplot(1,2,2);
q=semilogy(t,abs(err_rk2),t,abs(err_rk4));
xlabel('t');
ylabel(texlabel('absolut error'));
set(q(1),'Color','Blue','Linewidth',2);
set(q(2),'Color','Red','Linewidth',2);
legend(q,{'RK2','RK4'});
legend(q,'Location','East');
