Solving differential equations with Ode15s

Hello, I have the following system:

CSTR in gaseous phase, transient conditions and N components. For each component i the mass balance is:

dCi/dt=1/V (QeCi,e-QsCi,s )+γi r

where V is the volume of the reactor, Qe and Qs are the volumetric flowrates of the inlet and outlet respectively, Ci,e and Ci,s are the inlet and outlet concentrations of component i, and r is the rate of the reaction that takes place.

However Qs is unknown so I summed up all the equations above and got this:

Qs =Qe*(1/∑i Ci,s)*(∑i Ci,e -V/Qe ∑i dCi/dt -V/Qe ∑i γir)

I need to solve this set of equations, where Qe and Ci,e vary with time.

I defined the equations above in a function dydt=function(t,y,Ce,Qe,tspan) so I have N differential equations + equation Qs. I want to solve this set. I would also like to know Qs and not just the concentrations.

I use interp1 to have access to Ce and Qe values for each t. I defined the initial condition y0 where all values are zero except one of them (there is only one component entering the reactor).

I use [t,y]=ode15s(@(t,y)function(t,y,Ce,Qe,tspan),tspan,y0);

However I get a warning: Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.

I don't know why I get this error. Is my method correct and is it possible to have Qs as an output of the function?

[Torsten]

Qs can't be determined by adding up the equations, but only by an additional independent condition, e.g. that the volume remains constant or the pressure remains constant or something similar.

Best wishes

Torsten.

[FM]

The pressure and the temperature remain constant but I can't figure out what that independent condition would be. I use equatioin PQ=FRT so Qs=QeFs/Fe but Fs is also unknown. That's why I summed the equations, because I need an extra equation but I don´t know what it is.

I was trying to figure out where the problem comes from and I added a condition,

 if t==0
   Qs=0
 else 
   Qs= expression above
end

if I do that it no longer shows a warning.

So now I'm trying to understand how ode15s works. I thought the time t in the input of the dydt function corresponds to tspan given during ode15s no? Like if tspan is 0,1,2,3... 10 then t starts at 0 and then is equal to 1, 2 etc? But it seems as though it takes the value 0 multiple times before moving to the next...