Fixed-bed reactor problem

I want solve this system of non-linear ODEs.

rate equation is

.

and balances equation is

.

this is my function script

function f = ODEfun(Z,Y)
global bi
Xb = Y(1);
Xc = Y(2);
T = Y(3);

% PARAMETER
P = 1.1;
G = 4684;
hw = 134;
U = 82.7;
dHB = -307;
dHC = -1090;
ya0 = 0.0093;
yo0 = 0.208;
rhob = 1300;
Dp = 0.003;
D = 0.025;
cp = 0.25;
Mf = 29.48;
R = 8.31447;
Tj = 653.15;

k1 = exp((-27000/R*T)+19.837);
k2 = exp((-31000/R*T)+20.86);
k3 = exp((-28600/R*T)+18.97);

rb = (P^2)*yo0*ya0*(k1*(1-Xb-Xc)-k2*B);
rc = (P^2)*yo0*ya0*(k2*B+k3*(1-Xb-Xc));

dXbdZ = (Mf*rb*rhob)/(ya0*G*(1+bi));
dXcdZ = (Mf*rc*rhob)/(ya0*G*(1+bi));
dTdZ = (-4*U)/(D*G*cp)*(T-Tj)+(rhob*(rb/(1+bi))*(-dHB))/(G*cp)+(rhob*(rc/(1+bi))*(-dHC))/(G*cp);

f = [dXbdZ; dXcdZ; dTdZ];
end

clc
global bi
bi = 0;
Zspan = linspace(0, 3, 500);
y0 = [0 0 623];
[z y] = ode45(@ODEfun,Zspan,y0);

but result is strange.

Xb, Xc does not change.

Answer 1

[darova]

You made mistakes

I tried another Zspan

Zspan = [0 1e-6];
y0 = [0 0 623];
[z y] = ode45(@ODEfun,Zspan,y0);

And get this