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

%
%  We want to perform three simulations
%
hvec = [ 0.32 , 0.032 , 0.0032 ]

hvec =

    0.3200    0.0320    0.0032

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

h =

    0.0032

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);
[n,m] = size(ycorr);
nm = n*m;
%
%  We loop over the three values of h
%
%  Case (a):
%
h1 = hvec(1);
t1 = 0:h1:10;
x = x0;
y1 = c*x0;
[n,m] = size(t1);
n = n*m - 1;
theta = 0.45;
tol = 1.0e-5;
echo off
y = ycorr(1:100:nm)';
%
%  Compute the error
%
errabs = y - y1;
er1 = norm(errabs,'inf')

er1 =

    0.0912

errrel = errabs ./ y;
erl1 = norm(errrel,'inf')

erl1 =

    0.0021

%
%  Case (b):
%
h2 = hvec(2);
t2 = 0:h2:10;
x = x0;
y2 = c*x0;
[n,m] = size(t2);
n = n*m - 1;
theta = 0.45;
tol = 1.0e-5;
echo off
y = ycorr(1:10:nm)';
%
%  Compute the error
%
errabs = y - y2;
er2 = norm(errabs,'inf')

er2 =

  3.4248e-006

errrel = errabs ./ y;
erl2 = norm(errrel,'inf')

erl2 =

  1.5978e-007

%
%  Case (c):
%
h3 = hvec(3);
t3 = 0:h3:10;
x = x0;
y3 = c*x0;
[n,m] = size(t3);
n = n*m - 1;
theta = 0.45;
tol = 1.0e-5;
echo off
y = ycorr';
%
%  Compute the error
%
errabs = y - y3;
er3 = norm(errabs,'inf')

er3 =

  1.2194e-006

errrel = errabs ./ y;
erl3 = norm(errrel,'inf')

erl3 =

  1.6008e-007

%
%  Plot the results
%
subplot(3,1,1)
plot(t,ycorr,'k-',t1,y1,'r-.')
grid on
title('BI4/5(0.45) Integration with Different Step Sizes')
xlabel('Time')
ylabel('y')
subplot(3,1,2)
plot(t,ycorr,'k-',t2,y2,'r-.')
grid on
xlabel('Time')
ylabel('y')
subplot(3,1,3)
plot(t,ycorr,'k-',t3,y3,'r-.')
grid on
xlabel('Time')
ylabel('y')
print -deps nsds_hw3_14.eps
%
%  Close diary and exit
%
diary off