function [x,A,nrDoFs] = FEM10(c4n,n4e,n4sDb,f,u4Db)
    N=size(c4n,1); d=size(c4n,2);
    A=sparse(N,N); b=zeros(N,1); x=zeros(N,1);
    for j=1:size(n4e,1)
    	area=abs(det([ones(1,d+1);c4n(n4e(j,:),:)'])/factorial(d));
    	grads=[ones(1,d+1);c4n(n4e(j,:),:)']\[zeros(1,d);eye(d)];
    	A(n4e(j,:),n4e(j,:))=A(n4e(j,:),n4e(j,:))+area*(grads*grads');
    	b(n4e(j,:))=b(n4e(j,:))+f(sum(c4n(n4e(j,:),:))/(d+1))*ones(d+1,1)*area/(d+1); 
    end
    freeNodes=setdiff(1:N,unique(n4sDb)); nrDoFs=length(freeNodes);
    x(unique(n4sDb)) = u4Db(c4n(unique(n4sDb),:));
    b = b - A*x;
    x(freeNodes)=A(freeNodes,freeNodes)\b(freeNodes); 
end