function U=solveBVPsystem(A,b,ID,R,method)
  % ID: indices Dirichlet
  % R: relevement a 0 de gD
  % Usage:
  %   V=FEM.solveBVPsystem(A,b,[]);      % Pas de conditions de Dirichlet
  %   V=FEM.solveBVPsystem(A,b,ID,R);    % methode 2 par defaut
  %   V=FEM.solveBVPsystem(A,b,ID,R,1);  % methode 1
  if isempty(ID),U=A\b;return;end
  if nargin==4,method=2;end   
  assert(ismember(method,[1,2]))
  nq=length(b);
  switch method
  case 1
    I=speye(nq);
    A(ID,:)=I(ID,:);
    b(ID)=R(ID);
    U=A\b;
  case 2
    IDc=setdiff(1:nq,ID);
    U=R;
    U(IDc)=A(IDc,IDc)\(b(IDc)-A(IDc,ID)*R(ID));
  end
end
