I'm not a chemist, but are you sure your equations are correct? Looking at the change of B, it is initialized at zero, and then pulled down to minus infinity. The reason is that no matter the sign of y(1), B becomes smaller. And being initialized at zero, y(2) itself will always be smaller or equal to zero. There is no way this can be stable, hence the integration fails. Also, the difference in scale between the different terms is bound to bring numerical issues, even if the equations are correct.
On a sidenote: It's easier to read if you use e.g. dy and y than changing the names from y to A,B to F. Also, global variables are really not pretty in my opinion. A better option would be:
% script
clear all; clc
%butanotiol+heptano
%T=25ºC
Deff=0.046E-8 %m2/s
ro_b=167.9400;
L=4.27;
U=0.0024;
k=4.68E-2
Ca0=6.05e8;
%ODE
z=0:0.1:L;
z1=z';
y0=[Ca0; 0]
ode_react = @(z,y)react(z,y,U,Deff,ro_b,k);
[z,y]=ode15s(ode_react,z1,y0)
plot(z,y(:,1))
ylabel('Concentración Ca [mol/m3]')
xlabel('z [m]')
end
% and your function
function dy=react(z,y,U,Deff,ro_b,k)
dy = zeros(2,1);
dy(1)=y(2);
dy(2)=(U./Deff).*y(2)-k.*ro_b.*(y(1).^2)./Deff;
end
