%  This function solves homework problem [H2.1]
%  --------------------------------------------
%
echo on
diary nsds_hw2_1.dry
%
%  Start by entering the system
%
A = [ 1250 -25113 -60050 -42647 -23999
       500 -10068 -24057 -17092  -9613
       250  -5060 -12079  -8586  -4826
      -750  15101  36086  25637  14420
       250  -4963 -11896  -8438  -4756 ]
b = [  5 ;  2 ;  1 ; -3 ;  1 ]
c = [ -1 , 26 , 59 , 43 , 23 ]
d = 0
x0 = [ 1 ; -2 ;  3 ; -4 ;  5 ]
%
%  Let's check where the eigenvalues are
%
l = eig(A)
%
%  All eigenvalues are in the left half plane.
%  Determine step size that leads to marginal stability
%
hmarg = min( -2*real(l) ./ ( abs(l) .* abs(l) ) )
%
%  We want to perform five simulations
%
hvec = [ 0.1*hmarg , 0.95*hmarg , hmarg , 1.05*hmarg , 2*hmarg ]
%
%  Let us perform a simulation in Matlab, and use it as
%  the "correct" reference solution, relying on Matlab's
%  capability to solve this problem correctly.
%
h = hvec(1)
S = ss(A,b,c,d)
t = 0:h:10;
u = ones(size(t));
ycorr = lsim(S,u,t,x0);
%
%  We loop over the five values of h
%
echo off
%
%  Case (a):
%
h1 = h;
t1 = t;
x = x0;
y1 = c*x0;
[n,m] = size(t1);
n = n*m - 1;
for i=1:n,
   xdot = A*x + b;
   x = x + h1*xdot;
   y1 = [ y1 , c*x ];
end
%
%  Case (b):
%
h2 = hvec(2);
t2 = 0:h2:10;
x = x0;
y2 = c*x0;
[n,m] = size(t2);
n = n*m - 1;
for i=1:n,
   xdot = A*x + b;
   x = x + h2*xdot;
   y2 = [ y2 , c*x ];
end
%
%  Case (c):
%
h3 = hvec(3);
t3 = 0:h3:10;
x = x0;
y3 = c*x0;
[n,m] = size(t3);
n = n*m - 1;
for i=1:n,
   xdot = A*x + b;
   x = x + h3*xdot;
   y3 = [ y3 , c*x ];
end
%
%  Case (d):
%
h4 = hvec(4);
t4 = 0:h4:10;
x = x0;
y4 = c*x0;
[n,m] = size(t4);
n = n*m - 1;
for i=1:n,
   xdot = A*x + b;
   x = x + h4*xdot;
   y4 = [ y4 , c*x ];
end
%
%  Case (e):
%
h5 = hvec(5);
t5 = 0:h5:10;
x = x0;
y5 = c*x0;
[n,m] = size(t5);
n = n*m - 1;
for i=1:n,
   xdot = A*x + b;
   x = x + h5*xdot;
   y5 = [ y5 , c*x ];
end
%
%  Plot the results
%
echo on
title('Forward Euler Integration with Different Step Sizes')
subplot(3,2,1)
plot(t1,ycorr,'k-',t1,y1,'r-.')
grid on
xlabel('Time')
ylabel('y')
subplot(3,2,2)
plot(t1,ycorr,'k-',t2,y2,'r-.')
grid on
xlabel('Time')
ylabel('y')
subplot(3,2,3)
plot(t1,ycorr,'k-',t3,y3,'r-.')
grid on
xlabel('Time')
ylabel('y')
subplot(3,2,4)
plot(t1,ycorr,'k-',t4,y4,'r-.')
grid on
xlabel('Time')
ylabel('y')
subplot(3,2,5)
plot(t1,ycorr,'k-',t1,y1,t5,y5,'r-.')
grid on
xlabel('Time')
ylabel('y')
print nsds_hw2_1.eps -deps
%
%  Close diary and exit
%
diary off
return