%  This function solves homework problem [H3.15]
%  ---------------------------------------------
%
echo on
diary nsds_hw3_15.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 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
while t < 10,
  [x,err,t,h] = bi45tv_step(x,t,h,err,theta,tol);
  yvec = [ yvec , c*x ];
  tvec = [ tvec , t ];
end
echo on
[n,m] = size(tvec);
nm = n*m;
hvec = [ 0.1 , tvec(2:nm) - tvec(1:nm-1) ];
havg = sum(hvec) / (nm-1)
%
%  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
for i=1:nm,
  t = tvec(i);
  y = sum( r2 .* exp(l2*t) ) + c*expm(A*t)*x0;
  ycorr(i) = y;
end
echo on
%
%  Compute the errors
%
errabs = ycorr - yvec;
errrel = errabs ./ ycorr;
errmax = norm(errrel,'inf')
%
%  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
echo off
return