Problem with multigrid method implementation in MATLAB

I'm implementing a 2-grid method for solving discrete Poisson equation. I have read many papers to understand multigrid methods, for this implementation I'm following an example from this paper (p. 10). I believe I've implemented everything correctly, but I see a weird behaviour.

After prolonging the coarse grid correction back to fine grid the smooth error component should have become smaller, bringing us closer to the solution and at the same time the error should have become more oscillatory. I don't see this behaviour in my implementation. The error barely changes after coarse grid correction. I tried following another example (p. 38) and I see the same problem.

Am I not understanding something correctly? Here is my code and results:

%Grid size
m = 32 - 1;
h = 1/(m+1);
X = h:h:(1-h);

%The fine grid has m=31 points.
%The coarse grid has mr=15 points
mr = (m-1)/2;

%Restriction and prolongation matrices
%R - full weighting
R = zeros(mr, m);
for i = 0:mr-1
    R(i + 1, 2*i + 1) = 1/4;
    R(i + 1, 2*i + 2) = 1/2;
    R(i + 1, 2*i + 3) = 1/4;
end

%P - linear interpolation
P = zeros(m, mr);
for j = 0:mr-1
    P(2*j + 1, j + 1) = 1/2;
    P(2*j + 2, j + 1) = 1;
    P(2*j + 3, j + 1) = 1/2;
end

%------------------------------------------
%I'll be solving discrete Poisson equation
%-\nabla^2 u = f
%u(0)=u(1)=0

%Given function f
f = 25*pi^2*(sin(5*pi*X) + 9*sin(15*pi*X));

%System of equations Ax=b with initial guess x=0
%tridiag produces a square m by m matrix
A = tridiag(-1, 2, -1, m);
b = f' * h^2;
x = zeros(m,1);

%Analytic solution
u = sin(5*pi*X) + sin(15*pi*X);

%------------------------------------------
%Multigrid

%Pre-smoothing
%3 iterations of weighted Jacobi method with w=2/3
x = Jacobi(A,x,b,2/3,3);

%Up to this point everything is fine

%Restric the residue r -> rr
r = b - A * x;
rr = R * r;

%Ae = r on the coarse grid
Ar = tridiag(-1, 2, -1, mr);
er = zeros(mr,1);
%4 iterations of weighted Jacobi method with w=2/3
er = Jacobi(Ar,er,rr,2/3,4);

%Prolongation er -> e
e = P * er;
x = x + e;

function x = Jacobi(A, x0, b, w, Imax)
    n = size(x0, 1);
    x = zeros(n,1);
    I = 0;
    while I < Imax
        for i = 1:n
            x(i) = w/A(i,i) * (b(i) - sum(A(i,:)'.*x0) + A(i,i)*x0(i)) + (1-w)*x0(i);
        end
        x0 = x;
        I = I + 1;
    end
end