function MultDiffusB
clear, clc, format short g, format compact
tspan = [0 0.001]; % Range for the independent variable 
y0 = [0.0002229; 0; 0.007208]; % Initial values for the dependent variables function dYfuncvecdz = ODEfun(z,Yfuncvec); 
disp(' Variable values at the initial point ');
disp([' t    = ' num2str(tspan(1))]);
disp('           y                  dy/dt         ');
%disp([y0 ODEfun(tspan(1),y0)]);
%[t,y]=ode45(@ODEfun,tspan,y0);
NAB(:,1)=[2.396e-5; -3.363e-4];
disp([y0 ODEfun(tspan(1),y0,NAB(1,1),NAB(2,1))]);
err=1;
it=0;
while (err>1e-10) & (it<20)
    it=it+1;
    itno(it)=it;
    [t,y]=ode45(@ODEfun,tspan,y0,[],NAB(1,it),NAB(2,it));
    f(:,it)=[y(end,1); y(end,2)-2.701e-3];
    err=sqrt(f(:,it)'*f(:,it));
    for j=1:2
        delj=abs(NAB(j,it))*0.01;
        NAB(j,it)=NAB(j,it)+delj;
        [t,yp]=ode45(@ODEfun,tspan,y0,[],NAB(1,it),NAB(2,it));
        fp=[yp(end,1); yp(end,2)-2.701e-3];
        for k=1:2
            DF(k,j)=(fp(k)-f(k,it))/delj;
        end
        NAB(j,it)=NAB(j,it)-delj;
    end
    NAB(:,it+1)=NAB(:,it)-inv(DF)*f(:,it);
end
disp('           Iter. No.    NA     NB        f1            f2  ');
disp([itno' NAB(1,1:it)' NAB(2,1:it)' f(1,:)' f(2,:)']);
pause
for i=1:size(y,2)
    disp([' Solution for dependent variable y' int2str(i)]);
    disp(['              t                  y' int2str(i)]);
    disp([t y(:,i)]);
    plot(t,y(:,i));
    title([' Plot of dependent variable y' int2str(i)]);
    xlabel(' Independent variable (t)');
    ylabel([' Dependent variable y' int2str(i)]);
    pause
end
%- - - - - - - - - - - - - - - - - - - - - -
%function dYfuncvecdz = ODEfun(z,Yfuncvec);
function dYfuncvecdz = ODEfun(z,Yfuncvec,NA,NB); 
CA = Yfuncvec(1); 
CB = Yfuncvec(2); 
CC = Yfuncvec(3); 
%NB = -.0003363; 	%Molal flux of component B (kg-mol/m^2*s) 	
%NA = .00002396; 	%Molal flux of component A (kg-mol/m^2*s) 	
DAB = .000147; 	%Diffusivity of A through B (m^2/s) 	
NC = 0; 	%Molal flux of stagnant component A (kg-mol/m^2*s) 	
DAC = .0001075; 	%Diffusivity of A through C (m^2/s) 	
DBC = .0001245; 	%Diffusivity of B through C (m^2/s) 	
CT = .2 / (.082057 * 328); 	%Gas concentration g-mol/L 	
%errA = CA - 0; 	%Error in the calcaulated value of CA at the initial estimate (z = 0.001 m) 	
%errB = CB - .002701; 	%Error in the calcaulated value of CB at the initial estimate (z = 0.001 m) 	
xA = CA / CT; 	%Mole fraction of A 	
xB = CB / CT; 	%Mole fraction of B 	
xC = CC / CT; 	%Mole fraction of C 	
dCAdz = (xA * NB - (xB * NA)) / DAB + (xA * NC - (xC * NA)) / DAC; 	%Concentration of A (g-mol/L) 	
dCBdz = (xB * NA - (xA * NB)) / DAB + (xB * NC - (xC * NB)) / DBC; 	%Concentration of B (g-mol/L) 	
dCCdz = (xC * NA - (xA * NC)) / DAC + (xC * NB - (xB * NC)) / DBC; 	%Concentration of C (g-mol/L) 
dYfuncvecdz = [dCAdz; dCBdz; dCCdz];