Solve DAE with Ode15s

Now I'm trying to use Matlab ODE solver ode5i to solve a set of equations: M*dT+S*T=f; where M is singular and M, S and f remain constant in the call of ode15i in one overall simulation time-step (let's say my total simulation period is 10 days, and one time-step is maybe 20 minutes.).

The boundary condition is that every overall simulation time-step, after I call ode15i, I'll reassign initial condition for the solver and also the f will change.

Now I guess I'm not quite sure whether I'm calling the solver right or describe the equations right: When the documentation talks about Mass matrix, is it required that the elements inside are all 1 and 0s?

This is the function that I describe the equations: M*dT+S*T=f;

function dwdt=odefun(t,w)

global S ff M2

dwdt=-M2^-1*S*w+M2^-1*ff;

end

In the code above, I swept the original M matrix and filled 1 to original 0 places, then the M2 matrix has inverse. So the dwdt here should be equal to dT in the equations.

Then I call the ode to solver the T and tspan here is one overall simulation time-step:

for n=1:nstep-1,

    tspan = [0 timestep];

    ff=f(:,n);

    [Tt,Y] = ode15s(@odenotfun,tspan,y0);

    Y1=Y';

    T(:,n+1) = Y1(:,end);

    T(Nip(1),n+1)=Tamb(n+1);
    T(Nip(2),n+1)=Tsky(n+1);
    T(Nip(3),n+1)=Tg;

    y0 = T(:,n+1);

end

Question here, is it required to use the odeset here like this: options=odeset('Mass',M3); (For M3, I just assign every nonzero element to be 1 from the original M matrix).

The result that comes out of this is not quite reasonable for me, anyone knows where my approach might be wrong?

Thanks a lot!