关于张弛法编程的问题

离岸风

下面给的代码是我写的顺序张弛求解椭圆方程的代码，当超张弛系数（代码里c）为1的时候可以正常出结果（程序会运行10多秒），但是系数一旦比1大一点点，立马程序就会结束运行，而且出了结果也不对，我试过比1小的系数，也是可以出结果的，但是时间比系数是1的长，但是书上说这个系数在1~2之间啊，有没有大神帮忙解决下！！！

------------------------------------------------------------------------------------------------------------------------------------------------
program main
implicit none
!声明变量与常数
integer::irec,i,j,flag,a,b
integer,parameter::m=101,n=1001
real,parameter::pi=3.14159,d=pi/50.0,w=0.0001,c=1
real::g(m,n),r(m,n)
real::start,finish
call cpu_time(start)

!给上下边界赋值
do i=1,m
   g(i,1)=-sin((i-1)*d)
   g(i,n)=-sin((i-1)*d)*cos((n-1)*d)
end do

!给定侧边界值
do j=2,n
    g(1,j)=0
    g(101,j)=0
    end do

!给估计值
do j=2,n-1
    do i=1,m-1
        g(i,j)=-sin((i-1)*d)
    end do
end do

flag=1

!求helmholtz方程的近似解

        a=2
        b=2

        !利用同时张弛法进行迭代逼近
        do while(flag==1)
        
        !同时算第m次所有点的余差
        do j=2,n-1
            do i=2,m-1
                r(i,j)=g(i-1,j)+g(i+1,j)+g(i,j-1)+g(i,j+1)-(4+2*d*d)*g(i,j)-4*d*d*sin((i-1)*d)*cos((j-1)*d)
            end do
        end do
            !对所有点的余差r进行判断
             flag=0
            loop1:do j=2,n-1
                do i=2,m-1
                    if(abs(r(i,j)/(4+2*d*d))>w) then                        
                        flag=1
                        a=i
                        b=j      
                        exit loop1
                    end if
                end do               
            end do loop1

          !求算m+1次的估计值
            do j=b,n-1
                do i=a,m-1                  
                 g(i,j)=g(i,j)+c*r(i,j)/(4+2*d*d)  
                end do
            end do
            
        end do
          call cpu_time(finish)
write(*,*) finish-start
   !写.grd文件
   open(unit=11,file="G:\g(x,y)1.grd",form="unformatted",access="direct",status="unknown",recl=1)
        irec=1
   do j=1,n
            do i=1,m
                                write(11,rec=irec) g(i,j)
                        irec=irec+1
                end do
    end do
    close(11)

end program