function u=FEMSolveSym(nodes,gMat,gVec,n2d,gDFunc)

%[...] = FEMESolveSym (...)
%  solves the system of linear equations for a numerical solution of a PDE
%
%  u=FEMSolveSym(nodes,A,b,n2d,gDFunc)
%
% 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 in a symmetric, banded form (see SBSolve() )
% 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
%         it may also be given as a scalar value
% u   is the vector with the values of the solution
%
%see also ReadMesh, ShowMesh, FEMEquation, ShowSolution, FEMValue

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

ug=-SBSolve(gMat,gVec);
n=length(n2d);
u=zeros(n,1);


for k=1:n
  if n2d(k)>0
    u(k) = ug(n2d(k));
  else
    if is_scalar(gDFunc) u(k) = gDFunc;
       else              u(k) = feval(gDFunc,nodes(k,1:2));
    endif % scalar
  endif
endfor


endfunction
