function [Th,U,varargout]=solveP1(f,gD,N,varargin)
  % 3eme parametre en sortie optionnel: matrice de Masse
  % 4eme parametre en sortie optionnel: matrice de Rigidite
  p = inputParser;
  p.addParamValue('method', 2, @(x) ismember(x,1:2) ); 
  p.parse(varargin{:});
  meshfile=fc_oogmsh.gmsh.buildmesh2d('square4',N);
  Th=fc_simesh.siMesh(meshfile);

  % Calcul de Ah et b (sans conditions de Dirichlet)
  K=siMesh.Rigidite(Th);
  M=siMesh.Masse(Th);
  Ah=K;
  F=siMesh.eval(Th,f);
  b=M*F;

  % Prise en compte des conditions de Dirichlet (comme inconnues, methode 1)
  Dlab=1:4; % Dirichlet labels
  ID=siMesh.get_toGlobal(Th,'d',1,'labels',Dlab);
  R=siMesh.eval(Th,gD,'d',1,'labels',Dlab); % Relevement
  U=FEM.solveBVPsystem(Ah,b,ID,R,p.Results.method);
  if nargout>=3,varargout{1}=M;end
  if nargout==4,varargout{2}=K;end
end
