如何用有限差分方法计算格子波尔茨曼方程的导数？

我在做关于核反应堆堆芯的毕设，用的是格子波尔茨曼方法（LBM），在已经写出正确直角坐标系的模拟代码后，我尝试在该代码上直接修改为柱坐标LBM，二者区别在于柱坐标的源项有关于中子贮量率G的偏导项。我再尝试应用差分法算偏导之后，发现用tecplot所作出的图像是模糊不清晰的。

问题相关代码，请勿粘贴截图

代码如下

//计算一阶导数
for(j = 0;j < Nz;j++)
{

        for (i = 0; i < Nr; i++)
        {
        if(i == 0)
            Gr1[i][j] = (G[i + 1][j] - G[i][j]) / dr;
            else if (i == Nr - 1)
                Gr1[i][j] = (G[i][j] - G[i - 1][j]) / dr;
            else
                Gr1[i][j] = (G[i + 1][j] - G[i - 1][j]) / 2 * dr;
        }
    }

    //计算二阶导数
    
    for(j = 0; j < Nz; j++)
    {

        for (i = 0; i < Nr; i++)
        {
        if(i == 0)
            Gr2[i][j] = (G[i + 1][j] - G[i][j]) / dr * dr;
            else if (i == Nr - 1)
                Gr2[i][j] = (G[i][j] - G[i - 1][j]) / dr * dr;
            else
                Gr2[i][j] = (G[i + 1][j] + G[i - 1][j] - 2 * G[i][j]) /  dr * dr;
        }
    }

    for(i = 0; i < Nr; i++)
    {

        for (j = 0; j < Nz; j++)
        {
        if(j == 0)
            Gz2[i][j] = (G[i][j + 1] - G[i][j]) / dz * dz;
            else if (j == Nz - 1)
                Gz2[i][j] = (G[i][j] - G[i][j - 1]) / dz * dz;
            else
                Gz2[i][j] = (G[i + 1][j] + G[i - 1][j] - 2 * G[i][j]) / dz * dz;
        }
    }

运行结果及报错内容

所得数据在作图后是这样的

但我做出得直角坐标系得图是这样的

我的解答思路和尝试过的方法

由于我只增加了源项中的导数部分，所以我猜测是我导数没有写好

我想要达到的结果

想要得到的图像和图二一致

你的作图方式错了吧，图二是contour plot，而你的却是线条图