Weird behavior in pdepe at origin

I am have fluid dynamics model in spherical coordinates that I am simulating pdepe. It works well, except that I get weird bevahior at the origin, where it seems like pdepe is not maintaining zero-flux boundary conditions.

This figure shows the mass flux in the fluid near the origin at various times in the simulation:

You can see it is not zero at the origin. And if you look closely at the purple and green lines, you can see that if the mass flux was zero at the origin, the slope between r=0 and r=1 would more closely match the slope between r=1 and r=2.

The figure shows the mass density in the fluid near the origin:

These curves should have a zero derivative at the origin, but instead suddenly dip down near r=0. And, because of this, the mass density becomes negative at the origin at 40 ns (the purple curve), which is unphysical.

My questions are:

1) Is this expected/known behavior for pdepe?

2) Does the non-zero flux and/or negative density at the origin affect the solution? I know pdepe is attempting to set zero-flux boundary conditions at the origin, so does the value at the origin actually affect the rest of the solution?

3) Is there any way to fix this?

Thanks for your help!

[Torsten]

We neither know the equation(s) you try to implement nor the MATLAB code you are using. So we have no basis to answer your questions.

[Brian Cluggish]

Fair enough. The equations and code are pretty complicated. I'll see if I can reproduce the effect on a simpler model.

[William Rose]

@Brian Cluggish,

I agree that something is wrong when you have a negative concentration, and you can;t ignore it.

Maybe you can write the explicit forward-in-time difference equation, with centered spatial diferences on the right hand side. (See here.) You may have to use a quite small value for delta t to achieve stability. But it is probably the easiest difference equation to write and implement. You could at least test things out, even though it may be too ineficient for regular use.

When simulating diffusion-advection PDEs in cylindrical geometry, where v(r,t) the velocity of flow at radius r and time t, I use boundary conditions v(r,t)=0 when r=R (non-slip at the wall), and dv/dr=0 at r=0. And the radial flux=0 at r=0 and at r=R (at the center and at the wall).