ode15s convergence issue

Hello to everyone,

I am studying a transient reaction-diffusion system in the solid pellet with a detalied kinetics and I am having some (probably) numerical issues I d like to discuss with you.

My spatial variable is discretised using central finite difference method and my model function is structured as follows:

% calculation of reaction rates based on kinetic expressions (some of them non-linear i.e. k1*C1*C2^2-k2*C1
% calculation of concentrations time derivatives dcdt for all species in my
% model

However, I am getting somewhat disturbing (non physical) results in situation where kinetic rate constant k of two particular reactions I detected so far is set to some non-zero value whereas I am getting the concentration of species not involved in those two reactions different from zero (note that kinetic rate constants of reactions yielding conflicted species is set to zero, thus they cannot be produced physically)! Normally, I would doubt the accuracy of kinetic rate expressions in my code and/or reaction scheme but this seems not to be the case...

Concentrations are evaluated using ode15s considering that a system of equations to integrate is rather stiff due to kinetic rate expressions.

I did some tests and so far I detected some interesting behaviour:

1) ode45 though takes more time doesnt show numerical artefacts of conficted species.

2) variations of absolute and relative tolerances impact concentrations of conflicted species but those are still present with ode15s.

3) MaxStep , NonNegative, NormControl still yield conflicted results with ode15s

4) Reactions for which I observe those trends are non-linear in terms of kinetic expression

5) I am attaching those artefacts in this topic. Note that the order of magnitude of those artefacts is very small <10^-20.

6) Note that artefacts appear for larger time values.