% This program uses a generalized shooting technique in section 
% 7.2 to solve Eqs. (7.1.9), (7.1.10), (7.2.27), and (7.2.28)
% with boundary conditions (7.2.29). The guess of the missing 
% initial boundary condition is updated by Eq. (7.2.32).
%
max = 20;
CA0 = 0.9;
tol = 0.000001;

for I = 1:max,

    % initial conditions:
    x10 = 0; 
    dCA0 = 1;
    dx10 = 0;

    % solve the state and shooting systems simultaneously:
    c0 = [CA0 x10 dCA0 dx10]';
    clear t c;
    [t, c] = ode45n('model72',[ 0  1], c0);
    clear CA x1 dCA dx1;
    CA = c(:, 1);
    x1 = c(:, 2);
    dCA = c(:, 3);
    dx1 = c(:, 4);
    L = length(CA);

    % update CA0:
    delta = (1-CA(L)+.01*x1(L))/(-dCA(L)+.01*dx1(L));
    if abs(delta) < tol
	break;
    end
    CA0 = CA0 - delta;
end

iterations = I
t = 1 - t;
plot(t, CA);
title('Concentration Profile in the Turbular Reactor');
xlabel('Dimensionless Distance (z/L)');
ylabel('Demensionless Concentration (CA/CA0)');
%print fig76.eps;
