上述程序错误使用 \ 矩阵维度必须一致。

非常抱歉，程序中确实有一个错误。在使用 `\` 运算符时，矩阵维度必须一致。因此，我们需要将电位差矩阵 `U` 转换为列向量后再进行运算。下面是修正后的程序：

% 电阻率测深法理论曲线正演程序

% 输入参数
r = 50; % 电极距离
n = 20; % 电极数目
rho1 = 100; % 地下介质电阻率
rho2 = [500, 1000, 2000]; % 地下目标电阻率，多个值用逗号隔开
h = 2000; % 探测深度
dh = 100; % 深度间距

% 计算电极间距
d = r / (n - 1);

% 计算电极位置
elec_pos = (0:(n-1)) * d;

% 初始化电位数组
potential = zeros(n, n);

% 计算目标电位
for i = 1:n
    for j = 1:n
        if i ~= j
            distance = abs(elec_pos(i) - elec_pos(j));
            if distance <= r
                potential(i, j) = rho2(1);
            else
                potential(i, j) = rho1;
            end
        end
    end
end

% 有限差分法计算电阻率深度剖面
resistivity_profile = zeros(h/dh, numel(rho2));
for i = 1:numel(rho2)
    rho = rho2(i);
    for j = 1:(h/dh)
        z = j * dh;
        G = zeros(n-1, n-1);
        for k = 1:(n-1)
            G(k, k) = (rho1 + rho) / (2 * pi * (elec_pos(k+1) - elec_pos(k)) * z);
            for l = (k+1):n-1
                G(k, l) = -1 * (rho1 + rho) / (2 * pi * (elec_pos(l+1) - elec_pos(l)) * z);
                G(l, k) = -1 * (rho1 + rho) / (2 * pi * (elec_pos(l+1) - elec_pos(l)) * z);
            end
        end
        U = potential(2:end, 2:end) - potential(2:end, 1:end-1);
        I = G \ U(:);
        resistivity_profile(j, i) = (rho1 + rho) / (dh * sum(I));
    end
end

% 绘制电阻率深度剖面图
figure;
for i = 1:numel(rho2)
    loglog(resistivity_profile(:, i), (1:(h/dh)) * dh, 'LineWidth', 2);
    hold on;
end
xlabel('电阻率(Ωm)');
ylabel('深度(m)');
title('电阻率测深法理论曲线正演结果');
legend('rho2_1 = 500', 'rho2_2 = 1000', 'rho2_3 = 2000');

在这个程序中，我们将电位差矩阵 `U` 转换为列向量 `U(:)`，然后再进行运算，确保矩阵维度一致。