function CAFEM
c4n = load('c4n.dat'); n4e = load('n4e.dat'); n4Db = load('n4Db.dat');
dof4lvl = []; eta4lvl = []; energy4lvl = [];
while true
    [n4ed,ed4e,e4ed] = SIDES(n4e,n4Db);
    [x,A] = CFEM(c4n,n4e,n4Db);
    dof4lvl(end+1) = size(c4n,1);
    energy4lvl(end+1) = x'*A*x;
    eta4e = CESTIMATE(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 = x;
patch(X(n4e)',Y(n4e)',Z(n4e)',Z(n4e)');shading interp;
figure(2);
triplot(n4e,c4n(:,1),c4n(:,2));
figure(3);
loglog(dof4lvl,eta4lvl,'-x');hold all;
loglog(dof4lvl,sqrt( 0.2140758026865 - energy4lvl),'-o');
end