% This program uses the superposition principle to calculate the radial
% temperature gradients in an annular chemical reactor, described by
% Eq. (7.4.4) in section 7.4.
%

% generate one particular solution and two homogeneous solutions:
%
cp0 = [0 0]';
ch10 = [1 0]';
ch20 = [0 1]';
c0 = [cp0;
      ch10;
      ch20];
[r, c] = ode45n('model74',[ 0.03  0.04], c0);
Tp = c(:, 1);
DTp = c(:, 2);

plot(r, Tp);
axis([0.03 0.04 0 10]);
title('Particular Solution, Tp');
xlabel('Radius (ft)');
ylabel('Tp (degree F)');
print fig78.eps;
pause;

plot(r, DTp);
axis([0.03 0.04 0 2000]);
title('Particular Solution, DTp');
xlabel('Radius (ft)');
ylabel('DTp (degree F/ft)');
print fig781.eps;
pause;

Th1 = c(:, 3);
DTh1 = c(:, 4);
plot(r, Th1, '-', r, DTh1, '*');
axis([0.03 0.04 -0.1 1.2]);
title('First Homogeneous Solution')
xlabel('Radius (ft)');
ylabel('State Variables');
text(0.034, 1.1, 'Th1(r)');
text(0.034, 0.1, 'DTh1(r)');
print fig79.eps;
pause;

Th2 = c(:, 5);
DTh2 = c(:, 6);
plot(r, Th2);
axis([0.03 0.04 0.0 0.01]);
title('Second Homogeneous Solution')
xlabel('Radius (ft)');
ylabel('Th2(r)');
print fig791.eps;
pause;
plot(r, DTh2);
axis([0.03 0.04 0.0 1.0]);
title('Second Homogeneous Solution')
xlabel('Radius (ft)');
ylabel('DTh2(r)');
print fig792.eps;
pause;

% solve linear equation set to get the constants that satisfy the real
% boundary conditions.
%
LH1 = length(Th1);
LH2 = length(Th2);
Lp = length(Tp);

A = [Th1(LH1) Th2(LH2);
     DTh1(1)  DTh2(1)];

B = [800-Tp(Lp) 0-DTp(1)]';

const = A\B;

% compose the final solution. 
%
T = Tp + const(1)*Th1 + const(2)*Th2;

plot(r, T)
axis([0.03 0.04 790 800]);
title('Final Solution');
xlabel('Radius (ft)');
ylabel('T(r) (degree F)');
print fig710.eps;
