% Solution of second order ordinary differential equations
% Numerical method: modified Euler (Runge-Kutta, order 2)

clear all
close all

% ODE:
% a y''+ b y' + c y = 0
%
% Equivalent system of ODEs:
% v' = w
% w' = -b/a w - c/a v

a = 1;
b = 5;
c = 6;

Tmax = 10;
dt = 0.01;
t = 0:dt:Tmax;

v = zeros(1,length(t));
w = zeros(1,length(t));

v(1) = 2;
w(1) = 3;

for j=1:length(t)-1
    k1v = w(j);
    k1w = -b/a*w(j)-c/a*v(j);
    av = v(j)+k1v*dt;
    aw = w(j)+k1w*dt;
    k2v = aw;
    k2w = -b/a*aw-c/a*av;
    v(j+1)=v(j)+(k1v+k2v)*dt/2;
    w(j+1)=w(j)+(k1w+k2w)*dt/2;    
end

figure
hold on
plot(t,v,'-b','linewidth',2);
axis([0 Tmax -1 3]);
set(gca,'fontsize',20);
xlabel('t');
ylabel('y');




