%------------------------------------------------------------------------
%%%               ENGINEERING MATHEMATICS III - MATH 2Z03
%%%                               LAB #4
%%%                HIGHER ORDER ODEs: SPRING-MASS SYSTEM
%------------------------------------------------------------------------


%------------------------------------------------------------------------
% Covers:
% - "Advanced Engineering Mathematics" by D.G. Zill and M.R. Cullen,
%    Section 3.8.3 (example 6)
%------------------------------------------------------------------------


%------------------------------------------------------------------------
% Instructions:
% - Execute the commands given below in MATLAB. You can do this in
%   one of the following three ways:
%   (i) retype the commands one by one at the MATLAB command line (>>),
%   (ii) copy and paste the commands one by one to the MATLAB
%        command line (>>),
%   (iii) run all the commands at once by calling the present file
%         as a script, i.e. >> m2z03lab
% - Analyze the output and compare it with the comments embedded in
%   the file below.
%------------------------------------------------------------------------

clc; close all; clear all;

%------------------------------------------------------------------------
% Section 3.8.3: Spring-Mass Systems: Driven Motion
disp(' -----------------------------------------------------------')
disp(' Section 3.8.3: Spring-Mass Systems: Driven Motion ')

%%% Example:
% Solve the initial-value problem
%   0.2 x'' + 1.2 x' + 2x = 5cos(4t),

% Initial values:
%   1) x(0) = 0.5,   x'(0) = 0,
%   2) x(0) = 0.2,   x'(0) = 1,
%   3) x(0) = -1.5,  x'(0) = -4,
%   4) x(0) = 2,     x'(0) = 2.

% Exact solution (for the 1st initial conditions):
% x(t) = exp(-3t)*((38/51)*cos(t)-(86/51)*sin(t))-(25/102)*cos(4t)+(50/51)*sin(4t)

% After the substitution y' = u
% the equation is equivalent to the system:
%         x' = u,
%         u' = -6u - 10x + 25cos(4t)

% Heun_method_for_system

% Initial values
x01 = 0.5; u01 = 0;   % 1st initial value
x02 = 0.2; u02 = 1;   % 2nd
x03 = -1.5; u03 = -4; % 3rd
x04 = 2; u04 = 2;     % 4th

% Initial data
h = 0.05 % time step
T = 2    % maximum value of t
n = T/h  % number of steps

% Values of t, x and u at first step k = 1 (for different initial values)
t(1) = 0;
x1(1) = x01; u1(1) = u01; % 1st initial value
x2(1) = x02; u2(1) = u02; % 2nd
x3(1) = x03; u3(1) = u03; % 3rd
x4(1) = x04; u4(1) = u04; % 4th

% Loop to generate steps from k = 1 to k = n
for k = 1:n
  t(k+1) = t(k) + h; % increment for t
% 1st initial value
  x_star = x1(k) + h*u1(k);
  u_star = u1(k) + h*(-6*u1(k)-10*x1(k)+25*cos(4*t(k)));
  x1(k+1) = x1(k) + (h/2)*(u1(k)+u_star); % update for the solution x(t)
  u1(k+1) = u1(k) + (h/2)*(-6*u1(k)-10*x1(k)+25*cos(4*t(k))+ ...
			   -6*u_star-10*x_star+25*cos(4*t(k+1))); % update for the solution u(t)
% 2nd initial value
  x_star = x2(k) + h*u2(k);
  u_star = u2(k) + h*(-6*u2(k)-10*x2(k)+25*cos(4*t(k)));
  x2(k+1) = x2(k) + (h/2)*(u2(k)+u_star); % update for the solution x(t)
  u2(k+1) = u2(k) + (h/2)*(-6*u2(k)-10*x2(k)+25*cos(4*t(k))+ ...
			   -6*u_star-10*x_star+25*cos(4*t(k+1))); % update for the solution u(t)
% 3rd initial value
  x_star = x3(k) + h*u3(k);
  u_star = u3(k) + h*(-6*u3(k)-10*x3(k)+25*cos(4*t(k)));
  x3(k+1) = x3(k) + (h/2)*(u3(k)+u_star); % update for the solution x(t)
  u3(k+1) = u3(k) + (h/2)*(-6*u3(k)-10*x3(k)+25*cos(4*t(k))+ ...
         -6*u_star-10*x_star+25*cos(4*t(k+1))); % update for the solution u(t)
% 4th initial value
  x_star = x4(k) + h*u4(k);
  u_star = u4(k) + h*(-6*u4(k)-10*x4(k)+25*cos(4*t(k)));
  x4(k+1) = x4(k) + (h/2)*(u4(k)+u_star); % update for the solution x(t)
  u4(k+1) = u4(k) + (h/2)*(-6*u4(k)-10*x4(k)+25*cos(4*t(k))+ ...
         -6*u_star-10*x_star+25*cos(4*t(k+1))); % update for the solution u(t)
end

% Exact solution
xExact = exp(-3*t).*((38/51)*cos(t)-(86/51)*sin(t))-(25/102)*cos(4*t)+(50/51)*sin(4*t);

% Plotting (exact and numerical solutions)
plot(t,xExact,'-r',t,x1,'.-b'); % plotting (Figure 1)
title('Forced Damped Spring-Mass System: Solution'); % figure's title,
legend('exact solution','numerical approximation'); % legend
xlabel('t'); ylabel('x'); % and axis label

disp('  Press any key to continue ... ')
pause

% Plotting the phase portrait
figure; % opening new window for Figure 2
plot(x1,u1,'.-b'); hold on; % 1st trajectory
plot(x2,u2,'.-g'); hold on; % 2nd trajectory
plot(x3,u3,'.-r'); hold on; % 3rd trajectory
plot(x4,u4,'.-k'); hold off; % 4th trajectory
title('Forced Damped Spring-Mass System: Phase Portrait'); % figure's title,
legend(['x_{01} = ' num2str(x01) ', x\prime_{01} = ' num2str(u01)],...
       ['x_{02} = ' num2str(x02) ', x\prime_{02} = ' num2str(u02)],...
       ['x_{03} = ' num2str(x03) ', x\prime_{03} = ' num2str(u03)],...
       ['x_{04} = ' num2str(x04) ', x\prime_{04} = ' num2str(u04)]); % legend
xlabel('x(t)'); ylabel('x\prime(t)'); % and axis label

disp('  This is the end of the lab tutorial #4. ')