function u=FEMSolve(nodes,A,b,n2d,gD)

%[...] = FEMSolve (...)
%  solves the system of linear equations for a numerical solution of a PDE
%
%  u=FEMSolve(nodes,A,b,n2d,gD)
%
% nodes contains information about the mesh
%         see ReadMesh() for the description of the format
% A   is the matrix of the system to be solved.
%     It is stored as a full matrix
% b   is the RHS of the system to be solved.
% n2d is the renumbering of the nodes to the DOF of the system
%     n2d(k)=0  indicates that node k is a Dirichlet node
%     n2d(k)=nn indicates that the value of the solution at node k 
%               is given by u(nn)
% 'gD' is the function describing the Dirichlet boundary condition
%
% u   is the vector with the values of the solution
%
%see also ReadMesh, ShowMesh, FEMEquation, ShowSolution, FEMValue

if (nargin!=5) 
  help("FEMSolve");
  usage("FEMSolve(nodes,A,b,n2d,gD)");
endif

% solve the system of linear equations
ug=-A\b;
n=length(n2d);
u=zeros(n,1);

% evaluate the function on the Dirichlet section of the boundary and
% create a vector u with the solution
for k=1:n
  if n2d(k)>0
    u(k) = ug(n2d(k));
  else
    if is_scalar(gD) u(k) = gD;
       else          u(k) = feval(gD,nodes(k,1:2));
    endif % scalar
  endif
endfor
endfunction
