program differelye_laplace 
 use IMSL
 implicit none
 !差分方程的建立
    integer i,j,h,k,k1,k2,k3,k4,k5,k6,k7,k8
    real,parameter::lx=10,ly=10,lz=10
    real(kind=8) r(lx,ly,lz),c(lx,ly,lz),z(lx,ly,lz)
    real,parameter::pi=3.1415926
    integer,parameter::n=1000     !n=lx*ly*lz
    integer::dx,dy,dz,dx2,dy2,dz2
    real uu    
    real::a(n,n)=0
    real::b(n,1)=0
    real::u(n,1)=0

    dx=1/(lx-1)
    dy=2*pi/(ly-1)
    dz=1/(lz-1)
    dx2=dx*dx
    dy2=dy*dy
    dz2=dz*dz
     uu=2.0    

 do i=2,lx-1
  do j=2,ly-1
   do h=2,lz-1
    k=(h-1)*lx*ly+(i-1)*ly+j
       r(i,j,h)=dx*(i-1)
    c(i,j,h)=dy*(j-1)
    z(i,j,h)=dz*(h-1)
    a(k,k)=2/dz2+2/dx2+2/r(i,j,h)/dy2
    a(k,k-1)=1/dx2
    a(k,k+1)=1/dx2
    a(k,k-ly)=1/dx2*(1-dx/2/r(i,j,h))
    a(k,k+ly)=1/dx2*(1+dx/2/r(i,j,h))
    a(k,k-lx*ly)=1/r(i,j,h)/dy2
    a(k,k+lx*ly)=1/r(i,j,h)/dy2
   end do
  end do
 end do


 do i=1,lx
  do j=1,ly

    k1=(lz-1)*lx*ly+(i-1)*ly+j
    a(k1,k1)=1  
    b(k1,1)=10.0   
  end do
 end do

 do i=1,8
  do j=1,ly
  k2=(i-1)*ly+j
  a(k2,k2)=1
  b(k2,k2)=uu  
  end do
 end do

 do i=8,lx
  do j=1,ly
  k3=(i-1)*ly+j
  a(k3,k3)=1
  b(k3,k3)=u(lx*ly+k3,1)
  end do
 end do


 do j=1,ly
  do h=1,lz
   k4=(h-1)*lx*ly+(lx-1)*ly+j
   a(k4,k4)=1
   b(k1,1)=u((h-1)*lx*ly+(lx-2)*ly+j,1) 
   k8=(h-1)*lx*ly+j
   a(k8,k8)=1
   b(k8,1)=10*z(i,j,h)
  end do
 end do

 do i=1,lx
  do h=1,lz
   k5=(h-1)*lx*ly+(i-1)*ly+1
   k6=(h-1)*lx*ly+i*ly-1
   a(k5,k5)=1
   a(k6,k6)=1    
   b(k5,1)=b(k6,1)  
  end do
 end do

 do j=1,ly
  do h=1,lz
  k7=(h-1)*lx*ly+j
  a(k7,k7)=1
  b(k7,1)=b((h-1)*lx*ly+1,1)
  end do
 end do

  u=a.ix.b
  write(*,*) u
stop
end



总是提示：被零除，出现黄色箭头在第一个循环的a(k,k)=处