function [A,b,utilde] = FEM(P,E,BilinInts,BilinCoeffs,LinInts,LinCoeffs,ID,uhat,IN,tauhat)

d = size(P,1);
N = size(P,2); m = size(E,2); ND = length(ID); mN = length(IN); n = N-ND;
Atilde = sparse(N,N); btilde = zeros(N,1);
for e=1:m
    Ee = E(:,e); Pe = P(:,Ee);
    Ae = zeros(d+1,d+1); be = zeros(d+1,1);
    for k=1:length(BilinInts)
        Ae = Ae + feval(BilinInts{k},Pe,BilinCoeffs{k}(e));
    end
    for k=1:length(LinInts)
        be = be + feval(LinInts{k},Pe,LinCoeffs{k}(e));
    end
    Atilde(Ee,Ee) = Atilde(Ee,Ee) + Ae;
    btilde(Ee) = btilde(Ee) + be;
end
nonDirNodes = setdiff(1:N,ID);
A = Atilde(nonDirNodes,nonDirNodes);
utilde = zeros(N,1); utilde(ID) = uhat;
% btilde(IN) = btilde(IN) + tauhat;
btilde = btilde - Atilde*utilde; b = btilde(nonDirNodes);