function x = MFEM(c4n,n4e,ed4e,e4ed)
n = max(ed4e(:)); m = size(n4e,1);
B = zeros(3,3,m); C = zeros(3,m); b = zeros(m,1);
M = full(spdiags([ones(6,1),ones(6,1),2*ones(6,1),ones(6,1),ones(6,1)],...
    [-4,-2,0,2,4],6,6));
for j = 1 : m
    coord = c4n(n4e(j,:),:)';
    area = det([1,1,1;coord])/2;
    N = coord(:)*ones(1,3)-repmat(coord,3,1);
    signum = (-1).^(e4ed(ed4e(j,:),1)~=j);
    C(:,j) = signum.*[norm(N([5,6],2));norm(N([1,2],3));norm(N([1,2],2))];
    B(:,:,j) = C(:,j).*N'*M*N.*C(:,j)'/(48*area);
    b(j) =  -area;
end
J = ed4e(:,cumsum(ones(3,3)))'; K = ed4e(:,cumsum(ones(3,3))')';
B = sparse(J(:),K(:),B(:),n,n);
J = ed4e'; K = cumsum(ones(m,3))';
C = sparse(J(:),K(:),C(:),n,m);
x = [B,C;C',sparse(m,m)]\[zeros(n,1);b];
end