Solving boundary value problem such that we obtain nonnegative solutions

I treat a toy example to get my point across. In reality I have to deal with a much more complex model.

Let us consider a one dimensional boundary value problem using the BVP5c solver. Two liquids enter at different points, move in opposite directions and react with eachother. Liquid 1 enters at x = 0 and moves in the positive direction. Liquid 2 enters at x = 1 and moves in the negative direction. We are interested in the equilibrium distribution of the two concentrations over the domain [0,1].

In essence, the problem is that at some point (e.g. factor = 1e5, see Matlab code) the computed derivative makes the concentration of liquid 2 (moving in the negative direction) negative. This should not be possible. One could combat that by increasing the mesh size, however, that is not a solution that can be incorporated in my actual model due to the computational time.

How to make sure that both concentrations are nonnegative ?

____________________________

MATLAB CODE

____________________________

The script that runs the bvp5c solver

x = linspace(0,1,100);
yinit = [ 5 2 ];
factor = 1;

solinit = bvpinit(x,yinit);

sol = bvp5c(@(x,c)Concentrations(x,c,factor),@(ya,yb)BoundaryCond(ya,yb),solinit);

c = deval(sol,x);

____________________________

Boundary conditions

function res = BoundaryCond(ca,cb)

res = [ ca(1)-10; cb(2)-6 ];

____________________________

ODE functions

function dcdx = Concentrations(x,c,factor)

dcdx(1) = -(c(1) + c(2));
dcdx(2) = factor*(c(1) * c(2));

[Amy]

Hi Kay, did you ever manage to resolve this problem? I am facing the same difficulties at the moment.

Thanks,

Amy