echo on

%  This script will illustrate, using Matlab,
%  finding the zero input response of a system.
%
%  For example, consider the following differential
%  equation where f(t) is the input and y(t) is the
%  output (using the D-operator notation):
%
%    (D^5 + 2D^4 + 10D^3 + 7D^2 + D + 1)y(t) = f(t)
%
%  with initial conditions on the forth, third, second,
%  first, and zeroeth derivative being 0, 0, -1, 10, 
%  and 0, respectively.
%
% (Hit any key to continue)
pause
%
%  For the zero input response, enter characteristic
%  polynomial in matlab (i.e. vector of coefficients
%  on the D-operator in decending order)
%  
cep = [1 2 10 7 1 1];
%
% (Hit any key to continue)
pause
%  Find Roots (i.e. modes of the system)
rts = roots(cep)

%   These modes will be used later to determine
%   the form of the zero-state response
%
%   (Hit any key to continue)
pause
%
%   Set up initial condition vector for initial
%   conditions on the zero-th through 4-th order
%   derivatives.
%
ic = [0 10 -1 0 0];
%
%  Set up matrix elements for initial conditions
%  equation.  Based on the modes the form of the
%  zero-input solution can be written as
%
%  y(t) = c1*exp(rts(1)*t)+c1*exp(rts(2)*t)+ ... 
%           ... + c5*exp(rts(5)*t)
%
%  Note the exp(rts(k)*t) term can be complex
%  however, since y(t) is real, all imaginary
%  components should cancel out.  Also note
%  that the above form is only valid provided
%  there are no repeated roots.
%
% (Hit any key to continue)
pause

%  A system of equations can be set up to
%  solve for the constants (c1 ... c5) based
%  on the initial conditions (i.e. t=0).  The
%  result is:
%
%   0 = c1        + c2         +  ... c5
%  10 = c1*rts(1) + c2*rts(2)  + .....c5*rts(5)
%  -1 = c1*rts(1)^2 + c2*rts(2)^2  + ...c5*rts(5)^2
%   0 = c1*rts(1)^3 + c2*rts(2)^3  + ...c5*rts(5)^3
%   0 = c1*rts(1)^4 + c2*rts(2)^4  + ...c5*rts(5)^4
%
%   From which the matrix equation can be set up
%   as   ic = m*cof;  where ic is the column vector of
%   initial conditions, m is the matrix of modes, and cof is
%   the column vector of constants.  In matlab, cof can be 
%   solved by cof = inv(m)*ic or cof = m \ ic;
%
%   (Hit any key to continue)
pause
%  
%   Create m matrix:
%
m = [rts'.^0;rts'.^1;rts'.^2;rts'.^3;rts'.^4];
%
%  Solve System of Equations to get coefficients
%  the \ in matlab, is like multiplying the vector
%  ic by the inverse of m, but actually uses another
%  more stable method to solve the system of equations.
%  Note that ic is a row vector.  Must change it to a
%  column vector with the transpose operator '

cof = m \ ic'

%  (Hit any key to continue)
pause

%  To plot function, create a time axis
%  look at modes of system to get an idea of the
%  time interval to use and to decide over what
%  range to plot the data.
rts

%  Highest frequency is 2.9286/(2*pi) = about 1/3
%  So the sampling rate should be about 1/30
%  Also note that for the largest magnitude time
%  constant is about 1/0.0261 = about 50, so we
%  should plot it for about 3 to 5 time constants 200.
%
%  (Hit any key to continue)
pause
%
%  Create time axis

t = [0:1/30:200];

%  Evaluate y(t) at every time axis point

y = cof'*exp(rts*t);

%  And plot it:

figure(1)
plot(t,real(y))   %  Use the real() command to
                  %  eliminate imaginary 0's
xlabel('Seconds')
ylabel('Signal Amplitude')

%  Is this waveform what you would expect, based on
%  the system modes?

echo off