function [u,S,M]=P1FEM(T,f,uD,uN,c)

% function [u,S,M]=P1FEM(T,f,uD,uN,c)
%
% -\Delta u + c u = g     in Omega
%               u = uD    on Gamma_D
%    \partial_n u = uN.n  on Gamma_N
% Input:
%     T     : basic triangulation
%     f     : vectorized function of two variables (source)
%     uD    : vectorized function of two variables (Dirichlet B.C.)
%     uN    : vectorized function of two variables
%             with two components   (uN.n = Neumann B.C.)
%     c     : constant parameter
% Output:
%     u     : P1-FEM approximation
%     S,M   : stiffness and mass matrices

nC=size(T.coordinates,1);
nE=size(T.elements,1);
dir=unique(T.dirichlet);
free=(1:nC)'; free(dir)=[];
u=zeros(nC,1);
rhs=zeros(nC,1);

[S,M,b1,b2]=FEMmatrices(T,f,uN);
matrix = S+c*M;
u(dir) = uD(T.coordinates(dir,1),T.coordinates(dir,2));
rhs=b1+b2-matrix(:,dir)*u(dir);
u(free)= matrix(free,free)\rhs(free);

return