How can I use the finite difference method to solve a simple unconstrained optimal control problem?

Question

I want to solve this problem using finite difference method.

\begin{equation} \begin{cases} &\min_{y, u} \quad \mathcal{J} = \frac{1}{2}\left\lVert y-z \right\rVert^{2} + \frac{\alpha}{2}\left\lVert u \right\rVert^{2}\\ &\begin{array}{ccc} s.t.&-\Delta y = u\quad in\quad \Omega \\ &y = 0 \quad on\quad \partial \Omega\\ \end{array} \end{cases} \end{equation}

According to this paper, enter link description here it should have a second-order convergence. But my numerical experiment results show that u is second-order convergent, while y remains basically unchanged when I set $z = y^* = \sin(\pi x)\sin(\pi y)$. Here is my code. Is there any bug?

clc;
clear;
close all;

n = 32; 
h = 1/n;
alpha = 1e-8;

B = diag(repmat([4], 1, n-1))+diag(repmat([-1], 1, n-2), 1)...
                +diag(repmat([-1], 1, n-2), -1);
I = eye(n-1);
delta_h = cell(n-1);
for i = 1:n-1
   for j = 1:n-1
    if i==j
        delta_h(i,j)={-B};
    elseif i == j+1
        delta_h(i,j)={I};
    elseif i == j-1
        delta_h(i,j)={I};
    else 
        delta_h(i,j)={zeros(n-1)};
    end
   end
end
delta_h = cell2mat(delta_h)/h/h;

A = [-alpha*delta_h, -eye((n-1)^2); eye((n-1)^2), -delta_h];

f1 = zeros((n-1)^2, 1);
for i = 1:n-1
   for j = 1:n-1
      z_ex((i-1)*(n-1)+j) = fun_z(j*h, i*h);
   end
end
for i = 1:n-1
   for j = 1:n-1
      u_ex((i-1)*(n-1)+j) = fun_u(j*h, i*h);
   end
end

f2 = z_ex;
F = [f1' f2]';

res = A\F;
y_cal = res(1:(n-1)^2);
y_L2 = sqrt(sum((y_cal(:)-z_ex(:)).^2)*h*h);
u_cal = res((n-1)^2+1:end)/alpha;
u_L2 = sqrt(sum((u_cal(:)-u_ex(:)).^2)*h*h);

u_ex = reshape(u_ex, [n-1, n-1]);
u_cal = reshape(u_cal, [n-1, n-1]);
y_cal = reshape(y_cal, [n-1, n-1]);
y_ex = reshape(z_ex, [n-1, n-1]);

Besides, When I change the values of $z\ and\ y^*$, for example, $xy(1-x)(1-y)$, the error may not even converge.

Thanks for your answer.