function [x,A] = NCFEM(c4n,n4e,ed4e,ed4Db)
n = max(ed4e(:)); m = size(n4e,1);
A = zeros(3,3,m); b = zeros(n,1); x = zeros(n,1);
for j = 1 : m
    area = abs(det([1 1 1;c4n(n4e(j,:),:)'])/2);
    grads = [1 1 1;c4n(n4e(j,:),:)']\[zeros(1,2);eye(2)];
    A(:,:,j) = 4*area*(grads*grads');
    b(ed4e(j,:)) = b(ed4e(j,:))+area/3;
end
J = ed4e(:,cumsum(ones(3,3)))'; K = ed4e(:,cumsum(ones(3,3))')';
A = sparse(J(:),K(:),A(:),n,n);
dof = setdiff(1:n,ed4Db(:));
x(dof) = A(dof,dof)\b(dof);
end