function NCAFEM
c4n = load('c4n.dat'); n4e = load('n4e.dat'); n4Db = load('n4Db.dat');
dof4lvl = []; eta4lvl = [];
while true
    [n4ed,ed4e,e4ed,ed4Db] = SIDES(n4e,n4Db);
    x = NCFEM(c4n,n4e,ed4e,ed4Db);
    dof4lvl(end+1) = size(n4ed,1);
    eta4e = NCESTIMATE(c4n,n4e,n4ed,ed4e,e4ed,x);
    eta4lvl(end+1) = norm(eta4e);
    if dof4lvl(end) >= 1e4; break; end
    mark4e = MARK(eta4e,0.5);
    [c4n,n4e,n4Db] = BISECT(c4n,n4e,n4Db,n4ed,ed4e,mark4e);
end
figure(1);
X = c4n(:,1); Y = c4n(:,2); Z = (ones(3)-2*eye(3))*x(ed4e)';
patch(X(n4e)',Y(n4e)',Z,Z); shading interp;
figure(2);
triplot(n4e,c4n(:,1),c4n(:,2));
figure(3);
loglog(dof4lvl,eta4lvl,'-x');hold all;
end