% Modified version of the MATLAB template file SolvBVDif.m used to model
% the tubular reactor with axial dispersion in Example 33.1 of "A First 
% Course on Kinetics and Reaction Engineering."
%
function [ax_length,ax_conv] = Example_33_1_ad
    % Known quantities and constants (in consistent units)
    k = 1.2E-2; % /s
    CAfeed = 1.0e-3; % mol/m^3
    Peax = 2.0;
    us = 0.01; % m/s
    d = 0.004; % m
    L = 50*d;
    Dax = us*d/Peax;

    % Function that evaluates the derivatives
	function dydx = bvodes(x,y)
        CA = y(1);
        CA_prime = y(2);
        r = k*CA;
		dydx = [
			CA_prime
			(r + us*CA_prime)/Dax;
		];
    end % of internal function bvodes

	% Function that calculates the errors at the boundaries
	function res = bvs(y_at_start,y_at_end)
        CA_inlet = y_at_start(1);
        CA_outlet = y_at_end(1);
        CA_prime_inlet = y_at_start(2);
        CA_prime_outlet = y_at_end(2);
		res = [
			us*CA_inlet - Dax*CA_prime_inlet - us*CAfeed;
            CA_prime_outlet;
		];
    end % of internal function bvs

    % All the necessary internal functions are now defined

	% Set up the initial mesh
	x_range_low = 0.0;
	x_range_high = L;
	n_mesh_points = 100.;
	
	% The next line creates an array for the independent variable
	x = linspace(x_range_low,x_range_high,n_mesh_points);

	% Guesses
	yinit = [
		CAfeed
        k*CAfeed/us
	];

	% Create a structure containing the mesh and guesses
	solinit=bvpinit(x,yinit);

	% Solve the odes
	result = bvp4c(@bvodes,@bvs,solinit);
    
    ax_length = result.x;
    ax_conv = (CAfeed - result.y(1,:))*100./CAfeed;

end % of file Example_33_1_ad.m
