% Modified version of the MATLAB template file SolvBVDif.m used to model
% the tubular reactor with axial dispersion in Activity 33.1 of "A First 
% Course on Kinetics and Reaction Engineering."
%
function [z,concA] = Activity_33_1_ad(Dax,L)
    % Known quantities and constants (in consistent units)
    k = 0.05; % dm3/mol/min
    CA0 = 1.0; % mol/dm3
    CZ0 = 0.0; % mol/dm3
    us = 6.0; % dm/min

    % Function that evaluates the derivatives
	function dydx = bvodes(x,y)
        CA = y(1);
        CA_prime = y(2);
        r = k*CA^2;
		dydx = [
			CA_prime
			(us*CA_prime + r)/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_prime_inlet = y_at_start(2);
        CA_outlet = y_at_end(1);
        CA_prime_outlet = y_at_end(2);
		res = [
			us*CA_inlet - Dax*CA_prime_inlet - us*CA0;
            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 = [
		CA0
        k*CA0^2/us
	];

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

	% Solve the odes
	result = bvp4c(@bvodes,@bvs,solinit);
    
    % Return CA vs. z
    z = result.x;
    concA = result.y(1,:);

end % of file Activity_33_1_ad.m
