function varargout=solveBVPsystemV2(A,b,ID,R,varargin)
  % ID: indices Dirichlet
  % R: relevement a 0 de gD
  if isempty(ID)
    assert(nargout==1)
    varargout{1}=A\b;
    return
  end
  p = inputParser;
  p.addParamValue('method', 2, @(x) ismember(x,[1,2]));
  p.addParamValue('solvefun', @(A,b) mldivide(A,b) );
  p.parse(varargin{:});
  method=p.Results.method;
  solvefun=p.Results.solvefun;
  nq=length(b);
  varargout=cell(1,nargout);
  switch method
  case 1
    I=speye(nq);
    A(ID,:)=I(ID,:);
    b(ID)=R(ID);
    [varargout{:}]=solvefun(A,b);
  case 2
    IDc=setdiff(1:nq,ID);
    U=R;
    [varargout{:}]=solvefun(A(IDc,IDc) , b(IDc)-A(IDc,ID)*R(ID));
    U(IDc)=varargout{1};
    varargout{1}=U;
  end
end
