function Activity_38_1
%Activity_38_1 Code used in the solution of Activity 38.1 of AFCoKaRE

    % Known quantities and constants (in consistent units)
    L = 0.05; % cm
    k = 0.0015/60; % /s
    DeA = 1.8e-7; % cm^2/s
    CA = 1e-3; % mol/cm^3
    
    % Calculate analytically
    Thiele = L*sqrt(k/DeA)
    eta_analytic = tanh(Thiele)/Thiele

    % Define the internal functions that will be needed

	% Function that evaluates the derivatives
	function dydx = bvodes(x,y)
		dydx = [
			y(2)
            k/DeA*y(1)
		];
    end % of internal function bvodes

	% Function that calculates the errors at the boundaries
	function res = bvs(y_at_start,y_at_end)
        % y_at_start is a column vector containing y values at the starting
        % x boundary and y_at_end is a column vector containing y values at
        % the ending x boundary.
		res = [
			y_at_start(2)
            CA - y_at_end(1)
		];

    end % of internal function bvs

    % All the necessary internal functions are now defined

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

	% Guesses
	yinit = [
		CA/2
        CA/L
	];

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

	% Solve the odes
	result = bvp4c(@bvodes,@bvs,solinit);
    
    % Calculate any other desired quantities from the results
    npoints = length(result.x);
    eta_numeric = DeA*result.y(2,npoints)/L/k/CA

end % of file Activity_38_1.m
