diff --git a/examples/zzEstimatorWithScottZhang.m b/examples/zzEstimatorWithScottZhang.m
index 861cd38620697d82534a1273e1c10b63d5935fc9..98bdbd2cfc1f2ce21b9e9fc05e3a2502d799e6f2 100644
--- a/examples/zzEstimatorWithScottZhang.m
+++ b/examples/zzEstimatorWithScottZhang.m
@@ -10,7 +10,7 @@ pmax = 4;
 lamalg = 0.1;
 
 %% initialization for every polynomial degree
-[nElem, nDofs, errEst, iterations, time, directSolveTime] = deal(zeros(pmax, 1000));
+[nElem, nDofs, errEst, recovery, oscillation, iterations, time, directSolveTime] = deal(zeros(pmax, 1000));
 for p = 1:pmax
     %% setup geometry & spaces
     printLogMessage('*** p = %d (of %d) ***', p, pmax)
@@ -66,7 +66,7 @@ for p = 1:pmax
         while ~solverIsConverged
             solver.step();
             u.setFreeData(solver.x)
-            mu2 = estimateZZ(blf, lf, u);
+            [mu2, recoveryTerm, oscillationTerm] = estimateZZ(blf, lf, u);
             estimator = sqrt(sum(mu2));
             solver.tol = lamalg*estimator;
             solverIsConverged = solver.isConverged;
@@ -76,6 +76,8 @@ for p = 1:pmax
         %% estimate error and store data
         nDofs(p,ell) = getDofs(fes).nDofs;
         errEst(p,ell) = estimator;
+        recovery(p,ell) = sqrt(sum(recoveryTerm));
+        oscillation(p,ell) = sqrt(sum(oscillationTerm));
         iterations(p, ell) = solver.iterationCount();
         nElem(p,ell) = mesh.nElements;
         printLogMessage('number of dofs: %d, estimator: %.2e after %d iterations', ...
@@ -94,12 +96,21 @@ for p = 1:pmax
     end
 end
 
+%% plot
 figure()
+colors = [[0      0.4470 0.7410];
+          [0.8500 0.3250 0.0980];
+          [0.9290 0.6940 0.1250];
+          [0.4940 0.1840 0.5560];
+          [0.4660 0.6740 0.1880];
+          [0.3010 0.7450 0.9330];
+          [0.6350 0.0780 0.1840]];
 for p = 1:pmax
     idx = find(nDofs(p,:) > 0);
     complexity = [cumsum(nDofs(p, idx) .* iterations(p, idx))];
-    loglog(complexity, errEst(p,idx), '-o', 'LineWidth', 2, 'DisplayName', ['p=',num2str(p)]);
+    loglog(complexity, errEst(p,idx), '-o', 'LineWidth', 2, 'DisplayName', ['\mu for p=',num2str(p)], 'Color', colors(p,:));
     hold on
+    loglog(complexity, recovery(p,idx), '-x', 'LineWidth', 2, 'DisplayName', ['recovery for p=',num2str(p)], 'Color', colors(p,:));
 end
 for p = unique([1,pmax])
     x = complexity / complexity(1);
@@ -113,9 +124,11 @@ title('error estimator over number of total computational cost')
 
 %% local function for recovery based a posteriori error estimation
 
-function indicators = estimateZZ(~, ~, u)
+function [indicators, recoveryTerm, oscillationTerm] = estimateZZ(~, lf, u)
 
 p = u.fes.finiteElement.order;
+mesh =  u.fes.mesh;
+trafo = getAffineTransformation(mesh);
 qr = QuadratureRule.ofOrder(2*p);
 
 % Remark: Scott-Zhang operator for piecewise polynomials of degree p corresponds to nodal interpolation
@@ -130,7 +143,55 @@ gGradU2.setData(nodalInterpolation(CompositeFunction(@(grad) grad(2,:,:), Gradie
 
 % compute recovery term
 recoveryTermHandle = @(g1, g2, grad) (g1 - grad(1,:,:)).^2 + (g2 - grad(2,:,:)).^2;
-recoveryTerm = CompositeFunction(recoveryTermHandle, gGradU1, gGradU2, Gradient(u));
-indicators = integrateElement(recoveryTerm, qr);
+recoveryTermFunction = CompositeFunction(recoveryTermHandle, gGradU1, gGradU2, Gradient(u));
+recoveryTerm = integrateElement(recoveryTermFunction, qr);
+
+% element-wise integrals of residual and residual squared
+R = CompositeFunction(@(u,f) u(1,:,:) + u(4,:,:) + f, Hessian(u), lf.f);
+rIntegrals = integrateElement(R, qr);
+rSquaredIntegrals = integrateElement(CompositeFunction(@(R) R.^2, R), qr);
+
+% compute patches
+patches = computePatches(mesh.elements);
+patchAreas = mesh.computePatchAreas;
+
+% compute oscillation term (projection only for p=0, i.e., integral mean)
+oscillationTerm = zeros(1, mesh.nElements);
+boundaryVertices = mesh.edges(:,mesh.getCombinedBndEdges);
+interiorVertices = true(1, mesh.nCoordinates);
+interiorVertices(boundaryVertices) = false;
+interiorVertices = find(interiorVertices);
+for i = 1:length(interiorVertices)
+    z = interiorVertices(i);
+    rPatchIntegral = sum(rIntegrals(patches{z}));
+    rSquaredPatchIntegral = sum(rSquaredIntegrals(patches{z}));
+    diff_integral = rSquaredPatchIntegral - (rPatchIntegral)^2 / patchAreas(z);
+    oscillationTerm(patches{z}) = oscillationTerm(patches{z}) + diff_integral;
+end
 
+% combine terms suitably
+hT = sqrt(trafo.area);
+oscillationTerm = hT.^2 .* oscillationTerm;
+indicators = recoveryTerm + oscillationTerm;
 end
+
+%% local function for computing patches for each vertex
+
+function patches = computePatches(elements)
+    n = size(elements, 1);
+    N = size(elements, 2);
+
+    patchCounts = accumarray(elements(:), 1);
+    ptr = cumsum([1; patchCounts]);
+    patchElements = zeros(n*N, 1);
+
+    for i = 1:N
+        for k = 1:n
+            idx = elements(k,i);
+            patchElements(ptr(idx)) = i;
+            ptr(idx) = ptr(idx) + 1;
+        end
+    end
+
+    patches = mat2cell(patchElements, patchCounts);
+end
\ No newline at end of file