% filename Prob_7.m
clear all
clc
format short e
global L k Dab Cao alpha

L= 1e-3; Cao=0.2; k=1e-3; Dab=1.2e-9;
alpha=0;		% initial guess on shooting parameter
errr=1;
count=0;

%  use shooting method to solve ode.  Use Newton Raph to iterate on alpha
while errr>1e-12 &  count <100
	yo(1)=Cao; 
	yo(2)= alpha;
	yo(3)=0;
	yo(4)=1;
	tspan = [0 L];
    [z,y]=ode45('rhs7',tspan,yo);	% 4 order RK method from 0 to L 
%	For version 4 use
%	[z,y]=ode45('rhs7',0,L,yo);
	nn=size(y);
	n=max(nn);			% finds length of vector y
	if y(n,4)==0
	   'Newtons failed, derivative is zero'
 	    break
	end
	errr= y(n,2)/y(n,4);		% finds the values at endpoint of integration then
	alpha=alpha-y(n,2)/y(n,4);	%  adjust the shooting parameter alpha
	errr=abs(errr)
	count=count+1;
end

%  use analytical solution to compare to 4 order RK solution
for kk=1:n
	pos=z(kk);
	result(kk,1)=z(kk);
	result(kk,2)=anyl7(pos);
	result(kk,3)= y(kk,1);
	result(kk,4)= anyl7(pos)-y(kk,1);
end

' position   Analytical	 integrated   difference'
disp(result); 
count
y
plot(result(:,1),result(:,2),'r',result(:,1),result(:,3),'ob')
title('Diffusion/Reaction Inside Catalyst')
xlabel('x')
ylabel('c')