Dear Adam

It would of course allwayss be helpfull, from the decription of the system and the FEM problem you are trying to analyze. It would be more helpfull if I would know what kind of mesh generation you have used and over what kind/shape of domain.

The fact that the stiffness matrix is singular means of course that its rows/columns are linearly dependent and there is some problem with the basis functions and you might want to double check the code again and perhaps just try with different basis functions.

I am not completely familiar with Matlabs Fem procedures/algorithms, but if the stiffness matrix is singular then by definition it can not have an inverse. Sure if there is numerical instability Matlab can try to deal with it it using for example Singular Value Decomposition and removing the rows/columns '' corresponding to the singular value which are really close to zero. After you move the 'almost zero' singular values you will get a regular stiffness matrix which is best approximation of the original on such as its regular measured by Euclidean norm.

Also I am now referring to previous answer by Gregorz. If your system is an PDAE (PARTIAL DIFFERENTIAL ALGEBRAIC EQUATION) I am not sure if Matlab has readily designed algorithms for PDAE:s. And yes if you try to treat the PDAE eq. as ordinary PDE,  and form the stiffness matrix as you woud do for the PDE it will just be wrong because of the algebraic dependence of the variables that the ALGEBRAIC  eqs. indicate.

Again if I might sugest try to learn the 'prolongation procedure' or eaven better 'fast prolongation procedure' to modify your PDAE equations to INVOLUTIVE FORM. This is important !!! To put it simple: The procedure will find all the differential relations from algebraic ones to produce you a list of equations in ivolutive form (in object called Jet-Space) which can be solved used standard algorithms, which also exist in Matlab !

Of course if you would send me a more detailed description of the actual problem and equations I might give you a more detailed answer and offer better accistance !

If the mass matrix is singular, it means that your system of differential equations in fact is a system of differential algebraic equations, so called DAEs. In MATLAB you can solve it by means ode15s or ode23t solvers. Also you have to specify that your mass matrix is singular in odeset function:

http://www.mathworks.com/help/matlab/ref/odeset.html#f92-1023650

Hope this helps!

Dear Adam

It would of course allwayss be helpfull, from the decription of the system and the FEM problem you are trying to analyze. It would be more helpfull if I would know what kind of mesh generation you have used and over what kind/shape of domain.

The fact that the stiffness matrix is singular means of course that its rows/columns are linearly dependent and there is some problem with the basis functions and you might want to double check the code again and perhaps just try with different basis functions.

I am not completely familiar with Matlabs Fem procedures/algorithms, but if the stiffness matrix is singular then by definition it can not have an inverse. Sure if there is numerical instability Matlab can try to deal with it it using for example Singular Value Decomposition and removing the rows/columns '' corresponding to the singular value which are really close to zero. After you move the 'almost zero' singular values you will get a regular stiffness matrix which is best approximation of the original on such as its regular measured by Euclidean norm.

Also I am now referring to previous answer by Gregorz. If your system is an PDAE (PARTIAL DIFFERENTIAL ALGEBRAIC EQUATION) I am not sure if Matlab has readily designed algorithms for PDAE:s. And yes if you try to treat the PDAE eq. as ordinary PDE,  and form the stiffness matrix as you woud do for the PDE it will just be wrong because of the algebraic dependence of the variables that the ALGEBRAIC  eqs. indicate.

Again if I might sugest try to learn the 'prolongation procedure' or eaven better 'fast prolongation procedure' to modify your PDAE equations to INVOLUTIVE FORM. This is important !!! To put it simple: The procedure will find all the differential relations from algebraic ones to produce you a list of equations in ivolutive form (in object called Jet-Space) which can be solved used standard algorithms, which also exist in Matlab !

Of course if you would send me a more detailed description of the actual problem and equations I might give you a more detailed answer and offer better accistance !

I suppose you are working with DAE systems, as Grzegorz said you must use 0de 15, 

Suppose the following DAE system

x1dot=x1x2^2-x2

0=-x2+(3+x2x1^2)

You must define a mass

M=[1 0; 0 0];

Then if the DAE system is index 1 and 

%%[0 10] tspan; [0.4 0.1] initial conditions

options = odeset('Mass',M);

[t1,y1] = ode15s(@Exa3NLmodel,[0 10],[0.4 0.1],options);

plot(t1,y1)

%where @Exa3NLmodel is the nonlinear function  for example (different file)

function [xdot]=Exa3NLmodel(t,x)

xdot=[x(1)*(x(2)^2-1);

-x(2)+(3+x(2))*x(1)^2];

Note that the solver works only if the DAE system is index 1. Otherwise, you must transform the system of index n to index 1 if possible.  By considering the same procedure you can solve systems with n variables. 

Thanks guys for interaction.