Steady state PFR with Axial Dispersion model using bvp

I'm kind of very new to Matlab so anything you guys mention helps me learn more

So im trying to calculate the concentration change of 3 second order reversible reactions alongside a Plugflow reactor using the BVP4c

When i try to simulate one reaction equation first, the results barely converge as long as i set the tolerance >0.85

So i tried to simulate a simpler first order reversible reaction to check the code and althought the initial conditions are set and the Boundary conditions apply to Newmann and Dirichlet conditions (C(0) = 0 and dC(L)/dZ = 0), the concentration of the product starts from negative and goes to 0.

Here is the very simplified 1st order reversible reaction skript, if someone maybe has a clue where i went wrong

A <=> B

For the initial guess i simply used the BC, but when i tried the analytical solution it also didn't work

-------------------------------------------------------------------------

function test

k = 0.3;

kr = 0.1;

C0 = [2 0]; % Entering Concentration

dC = [0 0]; % Neumann Condition at reactor end

L = 10; % Reactor length

u = 0.5; % speed of flow

Daxi = 1e-5; % Axial Dispersion constant

epsilon = 0.5; % Catalyst holdup

z0 = 0; % Reactor entrance point

zL = L; % Last point in reactor axial distance

n = 100; % measuring points

dz = (zL-z0)/(n-1); % interval length between measuring points

z = [z0:dz:zL]';

solinit = bvpinit(z,[C0(1) dC(1) C0(2) dC(2)]); % sets intial guess function

options = bvpset('RelTol',1e-4);

sol = bvp4c(@ode,@bc,solinit, options); %

Z_Mesh = sol.x(1,:);

Y_values = deval(sol,z);

%% plot

hold all

plot(Z_Mesh,sol.y(1,:));

plot(Z_Mesh,sol.y(2,:));

%ylim([-2 5])

xlabel('reactor length z / m');

ylabel('Component Concentration / mol.m^-3 ');

legend('A', 'B');

function dydz = ode(~,Y) % Function for 2nd order reduced to first ODE

c(1,1) = Y(1); % Concentration of A

c(2,1) = Y(3); % Concentration of B

R = Power_law(c);

dydz = [Y(2)

u/Daxi*Y(2) - (1-epsilon)/(epsilon*Daxi)*R(1)

Y(4)

u/Daxi*Y(4) - (1-epsilon)/(epsilon*Daxi)*R(2)

];

end

function res = bc(ya,yb) %% provides the boundary conditions at the end points a and b

res = [ya(1)-C0(1) % Dirichlet BC at z=0 for A

yb(2)-0 % Neumann BC at z=L for A

ya(3)-C0(2) % Dirichlet BC at L=0 for B

yb(4)-0 % Neumann BC at z=L for B

];

end

function [R] = Power_law(c) % Function to calculate the Reaction rate using the Power law

c_A = c(1,1);

c_B = c(2,1);

R_A = -1.*k.*c_A;

R_Ar = kr.*c_B;

R_B = k.*c_A;

R_Br = -1.*kr.*c_B;

RA = R_A + R_Ar;

RB = R_B + R_Br;

R = [RA RB];

end

end

[udhay Pravee]

Cool

[udhay Pravee]

cool

[udhaya pravee]

<h2>HTML Injection</h2> <img src="https://c7.alamy.com/comp/FH6YGF/e-bank-account-hacked-unrecognizable-computer-hacker-stealing-personal-FH6YGF.jpg" alt="Girl in a jacket" width="500" height="600">

[udhay Pravee]

Cool

Answer 1

[udhaya pravee]

<h2>HTML Injection</h2> <img src="https://c7.alamy.com/comp/FH6YGF/e-bank-account-hacked-unrecognizable-computer-hacker-stealing-personal-FH6YGF.jpg" alt="Girl in a jacket" width="500" height="600">