%
%  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 us try an initial step size of h=0.1 (why not!)
%
t = 0;
tvec = t;
h = 0.1;
x = x0;
yvec = c*x0;
err = 0;
theta = 0.45;
tol = 1.0e-5;
echo off
[n,m] = size(tvec);
nm = n*m;
hvec = [ 0.1 , tvec(2:nm) - tvec(1:nm-1) ];
havg = sum(hvec) / (nm-1)

havg =

    0.3118

%
%  Get the analytical solution
%
S = ss(A,b,c,d);
G = tf(S);
Pu = 1;
Qu = [ 1 0 ];
U = tf(Pu,Qu);
Y = G*U;
[Py,Qy] = tfdata(Y,'v');
[r2,l2] = residue(Py,Qy);
ycorr = zeros(size(yvec));
echo off
%
%  Compute the errors
%
errabs = ycorr - yvec;
errrel = errabs ./ ycorr;
errmax = norm(errrel,'inf')

errmax =

  3.2851e-005

%
%  Plot the results
%
subplot(2,1,1)
plot(tvec,ycorr,tvec,yvec,'o')
grid on
title('BI4/5(0.45) Integration with Step Size Control')
xlabel('Time')
ylabel('y')
subplot(2,1,2)
stairs(tvec,hvec)
grid on
xlabel('Time')
ylabel('step size')
print -deps nsds_hw3_15.eps
%
%  Close diary and exit
%
diary off