% This program uses the quasilinearization method (which actually use
% superposition principle in a loop to approximate the real solution)
% to solve the nonlinear tubular reactor with dispersion in section 7.6
%
global Cn I;

% simulation parameters
%
max = 100;
tol = 0.0001;
Pe = 10;
R = 0.5;
C0 = 5.0;

Cn = C0 * ones(1,100);

hold off;

for J = 1:max,

    Cp0 = 0; DCp0 = 0; 
    Ch10 = 1; DCh10 = 0;
    Ch20 = 0; DCh20 = 1;

% use superposition principle to solve the problem one step forward
%
    for I = 1:100,
	c0 = [Cp0 DCp0 Ch10 DCh10 Ch20 DCh20]';
        clear z c
        [z, c] = ode45n('model76',[ 0  0.01], c0);
        L = length(z);
        Cp(I) = c(L, 1);
	Cp0 = Cp(I);
        DCp(I) = c(L, 2);
	DCp0 = DCp(I);
        Ch1(I) = c(L, 3);
	Ch10 = Ch1(I);
        DCh1(I) = c(L, 4);
	DCh10 = DCh1(I);
        Ch2(I) = c(L, 5);
	Ch20 = Ch2(I);
        DCh2(I) = c(L, 6);
	DCh20 = DCh2(I);
    end

    % solve linear equation set to get the constants that satisfy the real
    % boundary conditions.
    %
    A = [DCh1(1) DCh2(1);
         Ch1(100)+1/Pe*DCh1(100)  Ch2(100)+1/Pe*DCh2(100)];

    B = [0-DCp(1) C0-Cp(100)-1/Pe*DCp(100)]';

    const = A\B;

    % glue together the particular and homogeneous solutions to get the
    % next iteration results.
    %
    C = Cp + const(1)*Ch1 + const(2)*Ch2;
    DC = DCp + const(1)*DCh1 + const(2)*DCh2;

    zz = [0.01:0.01:1]';
    plot(zz, C, '-', zz, DC, '--');
    axis([0 1 0 7]);
    title('Convergence Sequence');
    xlabel('Dimensionless Distance (z)');
    ylabel('C(z) and DC(z)');
    text(0.2, 2.8, 'C1');
    text(0.1, 2.0, 'C2');
    text(0.1, 1.4, 'C3');
    text(0.6, 1.2, 'DC1');
    text(0.3, 1.2, 'DC2');
    text(0.02, 1.0, 'DC3');

    hold on;

    if sum(abs(C-Cn))<tol
	break;
    end
    Cn = C;

end
print fig712.eps
hold off
