Like this?:
global PT S DELH Q CPA CPB CPC CPD NA NC NB ND NE NT
PT = 1.5; % PT is total pressure in atmospheres
S = 0.05; % S= 0.5 square meters of area
DELH = 206014; % DELH is heat of reaction in Joules per gram mole
Q = 5000; % Heat input rate per unit lenth of bed (Joules/min-meter)
CPA = 36.9607; % specific heat (J/gmole-K)
CPB = 33.7295; % specific heat (J/gmole-K)
CPC = 29.1668; % specific heat (J/gmole-K)
CPD = 28.6455; % specific heat (J/gmole-K)
zz0 = 0.0;
zzf = 0.5;
YY0 = [3;3.5;0;0.000001;1300];
[zz,YY] = ode23(@ chrlspr3,[zz0 zzf],YY0);
K = 0.45; % For now, K is assumed to be independent of temperature, K = constant
NA= YY(:,1);
NB= YY(:,2);
NC= YY(:,3);
ND= YY(:,4);
T = YY(:,5);
NE= .45* YY(:,3).* YY(:,2)./YY(:,4);
NT= NA+NB+NC+ND+NE;
plot(zz,NA,'r+',zz,NB,'g-')
title('MOLES/MIN versus REACTOR LENGTH')
xlabel(' Reactor Length')
ylabel('Gram-moles/min')
grid
figure(2)
plot(zz,NC,'b-.',zz,ND,'y--',zz,NE,'r+')
title('MOLES/MIN versus REACTOR LENGTH')
xlabel(' Reactor Length')
ylabel('Gram-moles/min')
grid
figure(3)
plot(zz,T)
title('TEMPERATURE versus REACTOR LENGTH')
xlabel('Reactor Length')
ylabel('TEMPERATUTE 0K')
grid
function W = chrlspr3(~,YY)
global PT S DELH Q CPA CPB CPC CPD NT
S= 0.01;
NT = YY(1)+YY(2)+YY(3)+YY(4)+(0.45* YY(3)* (YY(2)/ YY(4)));
NTCP=YY(1)*CPA + YY(2)*CPB + YY(3)*CPC + YY(4)*CPD;
W3 = 1.02*(10^6)*S*(0.10949*PT*YY(1)/NT)/(1+0.239*YY(2)/YY(4)+17.62*YY(3)*PT/NT);
W =[-1.02*(10^6)*S*(0.10949*PT*YY(1)/NT)/(1+0.239*YY(2)/YY(4)+17.62*YY(3)*PT/NT);
-1.02*(10^6)*S*(0.10949*PT*YY(1)/NT)/(1+0.239*YY(2)/YY(4)+17.62*YY(3)*PT/NT);
W3
3*1.02*(10^6)*S*(0.10949*PT*YY(1)/NT)/(1+0.239*YY(2)/YY(4)+17.62*YY(3)*PT/NT);
(-S*DELH*W3+Q)/(NTCP)];
end
