Catalyst Modelling with MATLAB: Modeling Concentration Profile in Photoreactor

[AAMAIK]

I have little knowledge of MATLAB. I am trying to model the concentration profile along the photocatalyst film thickness in a photoreactor. Mass Transport takes place by diffusion and the rate of generation of component is given by rA.
Boundary conditions: Concentration of the component that is converted (phenol) on the surface of the catalyst is equal to the bulk
and the concentration gradient with respect to catalyst film thickness evaluated at the other end is equal to zero.

 

[jedishrfu]

Do you have a question?

 

[AAMAIK]

I am asking for links from which one could learn how to set up the boundary conditions and the differential equations in MATLAB to find a solution (concentration profile along the catalytic film). I have litttle experience in using MATLAB

 

[AAMAIK]

d^2c/dσ ^2+ θ^2*c/(1+Φ*c)=0
c evaluated at σ =0 is equal to 1
and the first derivative of c with respect to σ evaluated at σ =1 is equal to zero.
I put the equation above in a form of a system of first order equations by introducing the variables c1 and c2.
c1=c
c2=c'
I rewritten the governing equation as
c1'=c2
c2'=θ^2*c1/(1+Φ*c1)
with the transformed boundary condtions
c1(0)=1
c2(1)=0
The boundary value problem requires three pieces of information, the equation to be solved, its associated boundary condtions and initial guess for the solution and the parameters theta and phi.
phi=0.85
theta=4.5521
init=bvpinit(linspace(0,1,50), bc_init,phi,theta) % bc_init is an initial function profile
sol=bvp4c(@rhs_bvp,@bc_bvp, init);
σ= linspace(0, 1, 100)
BS=deval(sol,x)
plot(σ,BS(1, ;I want to build bc_init
function in=bc_init(σ);
in=[cos((pi/2)*σ)];
my input is going to be a vector is my guess for c1 and c2. Looking at the solutions at the boundary points I would guess the solution to behave like Cos(pi/2*σ)
So the solver will take my guess for the parameters phi and theta and produce an initial function evaluated at 50 points between 0 and 1 and send its on sol and operates on the equation and boundary and iterates towards a solution

function bc=bc_bvp(cl,cr,phi,theta)
bc=[cl(1)-1; cr(2)];
The boundary condition function of the vector c evaluated at the values σ=0,1

the function below represents the differential equation
function rhs=rhs_bvp(σ,c,phi,theta)
rhs=[c(2); ((theta)^2*c(1)/(1+phi*c(1)

I tried to run the program but its generating these errors
w = feval(v,x(1),extraArgs{:});

Error in Guess2 (line 5)
init2=bvpinit(linspace(0,1,50),@mat4init,phi,theta);