下面这个多重网格代码的误差太大，迭代次数也非常多，这说明它存在很多错误，请你按照多重网格的思想仔细检查是哪里出问题了，并修改我的代码，谢谢clc
clear
close all

%设置网格尺寸和大小
N = 64; 
L = 1;
h = L/N;
%u是一个N+1维的方阵，用来表示二维泊松方程的数值解
u = zeros(N+1, N+1);

% 计算真实解
x=0:h:L;
y=0:h:L;
[X, Y] = meshgrid(x, y);
u_true = sin(pi * X) .* sin(pi * Y); % 使用正确的元素乘法符，算出解析解

[X, Y] = meshgrid(x, y);
f = 2 * pi^2 * sin(pi *X).* sin(pi * Y);


error_iteration=zeros(1,1000);
for maxnum = 1:1000
    u = V_xunhuan(u,f,h);
    r = cancha(u,f,h);
    %计算误差并记录
    err=max(max(abs(u-u_true)));
    error_iteration(maxnum)=err;
    if err < 0.001
        break
    end
end
huatu_x=1:maxnum;
huatu_y=error_iteration(1:maxnum);
figure,plot(huatu_x, huatu_y,'r-');
title(['误差随迭代次数的变化 (一共迭代: ', num2str(maxnum), '次)']);
xlabel('迭代次数');
ylabel('误差');
grid on;

function jie = V_xunhuan(jie,f,h)
% 递归v循环
%预先平滑
jie = smoothing(jie,f,h);
%算残差
r_old = cancha(jie,f,h);
%限制
r_new = xianzhi(r_old);
%算误差
e = zeros(size(r_new));
%如果网格已经达到最粗，停止递归，否则继续递归
if length(e)-1 == 2
    e = smoothing(e,r_new,2*h);
else
    e = V_xunhuan(e,r_new,2*h);
end
% 修正
jie = jie + yantuo(e);
%后光滑
jie = smoothing(jie,f,h);
end


function res = smoothing(phi,f,h)
N = size(phi, 1) - 1;
res = phi;
for i = 2:N
    for j = 2:N
        res(i, j) = (phi(i+1, j)+res(i-1, j)+phi(i, j+1)+res(i, j-1)+h^2*f(i, j))/4;
      
    end
end
end

function res = cancha(phi,f,h)
N = size(phi, 1) - 1;
res = zeros(N+1, N+1);
res(2:N, 2:N) = f(2:N, 2:N)-(phi(1:N-1, 2:N)-2*phi(2:N, 2:N)+phi(3:N+1, 2:N)+phi(2:N, 1:N-1)-2*phi(2:N, 2:N)+phi(2:N, 3:N+1))/h^2;
end

function r_xi = xianzhi(r_cu)
N = (size(r_cu, 1) - 1) / 2;
r_xi = zeros(N+1, N+1);
for i = 2:N
    for j = 2:N   
       %直接映射r_xi(i, j) = r_cu(2*i, 2*j);
       %九点差分
      % r_xi(i, j) = 1/16 *( r_cu(2*i-1, 2*j-1)+r_cu(2*i-1, 2*j+1)+r_cu(2*i+1, 2*j-1)+r_cu(2*i+1, 2*j+1))+2/16*(r_cu(2*i, 2*j-1)+r_cu(2*i, 2*j+1)+r_cu(2*i-1, 2*j)+r_cu(2*i+1, 2*j))+4/16*r_cu(2*i,2*j);
       %五点差分 
       r_xi(i, j) = 2/16*(r_cu(2*i, 2*j-1)+r_cu(2*i, 2*j+1)+r_cu(2*i-1, 2*j)+r_cu(2*i+1, 2*j)) + 1/2 * r_cu(2*i, 2*j);      
    end
end
end


function e_xiwang = yantuo(e_cuwang)
    N = (size(e_cuwang, 1) - 1) * 2;  % 计算细网格的大小
    e_xiwang = zeros(N + 1, N + 1);  % 初始化细网格矩阵
    
    % 对角线（单元格中心）
    for i = 2:2:N
        for j = 2:2:N
            e_xiwang(i, j) =1/4*( e_cuwang((i+2) / 2, j / 2)+e_cuwang((i+2)/ 2, (j+2) / 2)+e_cuwang(i / 2, j / 2)+e_cuwang(i / 2, (j+2) / 2));
        end
    end
   
    % 原有点
    for i = 1:2:N + 1
        for j = 1:2:N + 1
            e_xiwang(i, j) = e_cuwang((i + 1) / 2, (j + 1) / 2);
        end
    end
    
    % 横向边中点
    for i = 1:2:N+1
        for j = 2:2:N
            e_xiwang(i, j) =1/2* (e_cuwang((i + 1) / 2, j / 2) + e_cuwang((i + 1) / 2,( j / 2) + 1));
        end
    end
    
    % 纵向边中点
    for i = 2:2:N
        for j = 1:2:N+1
            e_xiwang(i, j) = 1/2*(e_cuwang(i / 2, (j + 1) / 2) + e_cuwang((i / 2) + 1, (j + 1) / 2));
        end
    end
end