function [ gradDiff, toRefine ]  = z2indicator( mesh, u, epsTol, maxNum)

[ graduTri, gradu ]  = gradientReconstruction( mesh, u );

gradDiff = zeros(mesh.NuTri,1);
toRefine = zeros(mesh.NuTri,1);

for i=1:mesh.NuTri
  index    = mesh.T(i,:);
  errorX = sum(gradu(1,index) - graduTri(1,i))/3/2*mesh.Bdet(i);
  errorY = sum(gradu(1,index) - graduTri(1,i))/3/2*mesh.Bdet(i);
  gradDiff(i) = sqrt(errorX^2 + errorY^2);
end

[~, index] = sort(gradDiff);
toRefine(index(length(index)-maxNum:length(index)))=1;
toRefine(gradDiff<epsTol)=0;
