% SOLVEEXAMPLE--Example of mechanical system solved using the 
% symbolic toolbox of Matlab.

echo on % to display the progress of the program in the command window

% Define all variables and parameters that will appear in equations

syms Z1 Z2 F z2_0
syms M k s

% The equations will be written in a variable 'e'. This program is most
% easily modified if any previous version of 'e' is deleted from Matlab's
% workspace.

clear('e'); % removes 'e' from workspace if previously defined

% These are the equations of the system in the Laplace domain:

e{1} = F == k*(Z1-Z2);
e{2} = k*(Z1-Z2) == M*s^2*Z2 - M*s*z2_0 + k*Z2;

u = {Z1, Z2}; % these are the unknowns

% It is assumed that z2_0, denoting z2(0), and F are known. 
% The constants k and M are also assumed to be known.

% The following solves the system of equations for Z1 and Z2

x = [e, u];
a = solve(x{:});

% This instruction prints the calculated Z1

pretty(a.Z1)

% Another way to print Z1

a.Z1

% If f = f0, we can substitute F by f0/s:

syms f0
expr = subs(a.Z1, F, f0/s);

expr = simplify(expr) % rewrite expr to a form easier to read

% Let's tell Matlab that k and M are positive and real; by default they
% are assumed to be complex

syms k M positive

% Now, let's apply the inverse Laplace Transform

z1t = ilaplace(expr)

% Let's now assume some numerical values: M = 1kg, k = 55 N/m, z2_0 = 
% 5E-3 m and f0 = 2 N. The command below substitutes the numerical values:

z1t = subs(z1t, {M, k, z2_0, f0}, [1,55,5E-3,2]);

% Generally, z1t may not be a simple enough expression. It could be
% simplified in several steps, such as the following:

z1tn = vpa(z1t); % evaluate numerically all coefficients
% Uncomment the following line if using an old version of the symbolic toolbox (2013 or earlier)
% [z1tn, h] = simple(z1tn); 
z1tn = vpa(simplify(z1tn),3) % display the result with 3 significant digits

% Now we can plot the response of the system

t = 0:0.01:6; % time between 0 and 6 with interval 0.01 between points
z1 = subs(z1tn,'t',t);  % finds the numerical values of z1(t)
% You may have trouble with the line above if you use an
% old version of Matlab. If you have Maple, you could use the following 
% line instead:
% z1 = subs(maple('evalf',z1t),'t',t);
plot(t,z1) % plots the graph
title('z_1(t) for t > 0.'); % writes the title of the graph
ylabel('z_1(t) [m]'); % lables the ordinate
xlabel('t [s]'); % labels the abscissa

echo off

% He is wise in heart, and mighty in strength: who hath hardened himself 
% against him, and hath prospered? Job 9:4
