大家好，我最近在解一个2维的非线性扩散方程（主要是源项非线性），解法中用到ADI算法，每一个独立的方向又用追赶法解3对角矩阵。但解出来的结果非常失真，出现负值，并且越来越大，不知道是什么原因。希望大家伙出点建议，可能你们的一个点子就解决了。
下面是 追赶法在X方向的求解程序，大家帮忙看看：
         do z=2,Zb
       f(2,z)=a3(2,z)/a2(2,z)
       fe(2,z)=a2(2,z)
        do x=3,N-2
           fe(x,z)=a2(x,z)-(a1(x-1,z)*f(x-1,z))
           f(x,z)=a3(x,z)/fe(x,z)
        end do
        fe(N-1,z)=a2(N-1,z)-(a1(N-1,z)*f(N-2,z))
!----------------------------------------------------------
        y(2,z)=(b1(2,z)-a1(2,z)*Ud(1,z))/a2(2,z)                
        do x=2,N-1
           y(x,z)=(b1(x,z)-a1(x-1,z)*y(x-1,z))/fe(x,z)
        end do
        y(N-1,z)=(b1(N-1,z)-a3(N-1,z)*Sd(N,z)-a1(N-1,z)*y(N-2,z))/fe(N-1,z)
!----------------------------------------------------------
        Ud(N-1,z)=y(N-1,z)
    
        do x=N-2,2,-1
           Ud(x,z)=y(x,z)-f(x,z)*Ud(x+1,z)
        end do
!-----------------------------------------------------------
        Ud(N,z)=(4.0d0*Ud(x,z)+4.0d0*Ud(N-1,z)-Ud(N-2,z)-Ud(3,z))/6.0d0
        Ud(1,z)=Ud(N,z)
    end do
    
    do x=1,N
       Ud(x,1)=(4.0d0*Ud(x,2)-Ud(x,3))/3.0d0
       Ud(1,1)=Ud(N,1)
    end do