function MultDiffusA
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);
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);
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 
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];