%
%  Start by entering the system
%
A = [ 1250 -25113 -60050 -42647 -23999

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

       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 ]

b =

     5
     2
     1
    -3
     1

c = [ -1 , 26 , 59 , 43 , 23 ]

c =

    -1    26    59    43    23

d = 0

d =

     0

x0 = [ 1 ; -2 ;  3 ; -4 ;  5 ]

x0 =

     1
    -2
     3
    -4
     5

%
%  Let's check where the eigenvalues are
%
l = eig(A)

l =

  -4.0000 + 3.0000i
  -4.0000 - 3.0000i
  -5.0000          
  -2.0000          
  -1.0000          

%
%  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) ) )

hmarg =

    0.3200

%
%  We want to perform five simulations
%
hvec = [ 0.1*hmarg , 0.95*hmarg , hmarg , 1.05*hmarg , 2*hmarg ]

hvec =

    0.0320    0.3040    0.3200    0.3360    0.6400

%
%  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)

h =

    0.0320

S = ss(A,b,c,d)
 
a = 
                x1           x2           x3           x4           x5
   x1         1250  -2.511e+004  -6.005e+004  -4.265e+004    -2.4e+004
   x2          500  -1.007e+004  -2.406e+004  -1.709e+004        -9613
   x3          250        -5060  -1.208e+004        -8586        -4826
   x4         -750    1.51e+004   3.609e+004   2.564e+004   1.442e+004
   x5          250        -4963   -1.19e+004        -8438        -4756
 
b = 
       u1
   x1   5
   x2   2
   x3   1
   x4  -3
   x5   1
 
c = 
       x1  x2  x3  x4  x5
   y1  -1  26  59  43  23
 
d = 
       u1
   y1   0
 
Continuous-time model.
t = 0:h:10;
u = ones(size(t));
ycorr = lsim(S,u,t,x0);
%
%  We loop over the five values of h
%
echo off
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