% Example 12.11 system of equations (2nd order ODE)

clear all
close all
clc

% Set the initial parameters
x0 = [-10; 0];
h = 0.1;
t = [0:h:20];

% Exact result
p = 1;
zeta = 0.15/2;
wd = sqrt(1-zeta^2)*p;

xexact = @(t) x0(1).*exp(-zeta*p.*t).*(cos(wd.*t)+(zeta/wd).*sin(wd.*t));
vexact = @(t) -x0(1).*exp(-zeta*p.*t).*(p/sqrt(1-zeta^2)).*sin(wd.*t);

% Define the function
f =@(x,t) [x(2); -0.15.*x(2) - x(1)];

% Using explicit euler method
[tout,xout] = explicit_euler2(@(x,t) f(x,t),t,x0);


% Plot Result
figure
hold on
fplot(xexact,[min(t),max(t)],'LineWidth',2)
plot(tout,xout(:,1),'ok','MarkerFaceColor',[0,0,0],'MarkerSize',3)
plot(xlim,[0,0],'-k')
xlabel('t')
ylabel('x')
title(['Position: Explicit Euler with h = ',num2str(h)])
hold off

figure
hold on
fplot(vexact,[min(t),max(t)],'LineWidth',2)
plot(tout,xout(:,2),'ok','MarkerFaceColor',[0,0,0],'MarkerSize',3)
plot(xlim,[0,0],'-k')
xlabel('t')
ylabel('dx/dt')
title(['Velocity: Explicit Euler with h = ',num2str(h)])
hold off


% Using Runge-Kutta method
[tout2,xout2] = rk4sys(@(x,t) f(x,t),t,x0);


% Plot Result
figure
hold on
fplot(xexact,[min(t),max(t)],'LineWidth',2)
plot(tout2,xout2(:,1),'ok','MarkerFaceColor',[0,0,0],'MarkerSize',3)
plot(xlim,[0,0],'-k')
xlabel('t')
ylabel('x')
title(['Position: Runge-Kutta with h = ',num2str(h)])
hold off

figure
hold on
fplot(vexact,[min(t),max(t)],'LineWidth',2)
plot(tout2,xout2(:,2),'ok','MarkerFaceColor',[0,0,0],'MarkerSize',3)
plot(xlim,[0,0],'-k')
xlabel('t')
ylabel('dx/dt')
title(['Velocity: Runge-Kutta with h = ',num2str(h)])
hold off

