How to solve the 'DAE appears to be of index greater than 1' error for this simple system of equations?

Hey guys,

I am trying to solve a set of equations using ode15s and a mass matrix. However, I keep getting the error:

Error using daeic12 (line 76)
This DAE appears to be of index greater than 1.

I know the problem is in the mass matrix, when I turn it off the error dissapears. However, I have no idea how to solve this anymore, increasing the tolerances also has no influence.

The used code:

%% Call to ODE solver

par.numC = 3;

% Set integration to the radius of the catalytic particle

par.zspan = [0 par.rp];

% Create mass matrix: 0 for solid phase balances; 1 for fluid phase

M = spdiags([zeros(par.numC,1);ones(par.numC,1)],0,(2)*par.numC, (2)*par.numC);

% ODE settings

options=odeset('Mass',M,'RelTol',1e-6,'AbsTol',1e-6);

% options=odeset('RelTol',1e-6,'AbsTol',1e-6);

par.k_0 = 6.7404e-04; % Rate Constant [1/s]

par.cA_0 = 100; % Initial Concentration of A [mol/m3]

par.cB_0 = 100; % Initial Concentration of B [mol/m3]

par.cC_0 = 1e-9; % Initial Concentration of product C [mol/m3]

par.InitCond = repmat([par.cA_0; par.cB_0; par.cC_0;],2,1);

% Call to solver bv

[z,u] = ode15s(@EqTest,par.zspan,par.InitCond,options,par);

%% Plotter

par.r = linspace(0, par.rp, 2); % Create the cell face positions again

Conc_A = u(:,1:par.numC:end-2);

Conc_B = u(:,2:par.numC:end-1);

Conc_C = u(:,3:par.numC:end);

plot(par.r,Conc_A(end,:),par.r,Conc_B(end,:), par.r,Conc_C(end,:))

legend('A','B','C')

And the equations:

function [dudz,par]= EqTest(~,u,par)

% Reshape to column vector for ODE15s

u = reshape(u, par.numC, 2);

% Create a storage

dudz = zeros(size(u));

% Determine the kinetic rate

par.R1 = par.k_0 .* u(1,1:end-1) .* u(2,1:end-1) ;

% Calculate the concentration profile

% Component A

dudz(1,1:end-1) = - par.R1;

% Component B

dudz(2,1:end-1) = - par.R1;

% Component C

dudz(3,1:end-1) = + par.R1;

dudz = dudz(:) ;

end

Can anybody help me with this?

Thanks in advance