function DGAFEM
c4n = load('c4n.dat'); n4e = load('n4e.dat'); n4Db = load('n4Db.dat');
dof4lvl = []; eta4lvl = []; energy4lvl = [];
eta = 12; %penalty parameter
while true
    [n4ed,ed4e,e4ed,ed4Db] = SIDES(n4e,n4Db);
    [x,A,dof4e] = DGFEM(c4n,n4e,n4ed,e4ed,eta);
    dof4lvl(end+1) = numel(x);
    energy4lvl(end+1) = x'*A*x;
    eta4e = DGESTIMATE(c4n,n4e,n4ed,ed4e,ed4Db,e4ed,x,dof4e,eta);
    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(dof4e)',Z(dof4e)'); shading interp;
figure(2);
triplot(n4e,c4n(:,1),c4n(:,2));
figure(3);
loglog(dof4lvl,eta4lvl,'-x');hold all;
loglog(dof4lvl,sqrt( abs(0.2140758026865 - energy4lvl)),'-o');
end