cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  laplace.f: Solution of Laplace equation with finite differences     c
c								       c
c  comment: Output data is saved in 3D grid format used by gnuplot     c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      	Program laplace
      	Implicit none
      	Integer maxx, maxy
      	Parameter(maxx=40,maxy=40)
      	Real*8 x, p(maxx,maxy), pnew, dp, dpmax, omega
      	Integer i, j, iter, y
c open output file
      	Open(8, File='laplace.dat', Status='Unknown')
      	Open(9, File='laplace_dp.dat', Status='Unknown')

   	Do 5 i=1, maxx
          Do 6 j=1, maxy
   	   p(i,j)=0.0
 6       Continue   
 5     Continue   

c the side with constant potential   	
   	Do 10 i=1, maxx
   	   p(i,1)=100.0
 10 	Continue   
c
c iteration algorithm
  	Do 20 iter=1, 1000
          dpmax=0.0
  	  Do 30 i=2,(maxx-1)
  	    Do 40 j=2,(maxy-1)
  	      pnew=0.25*(p(i+1,j)+p(i-1,j)+p(i,j+1)+p(i,j-1))
              dp=pnew-p(i,j)
              dpmax=max(abs(dp),dpmax)
  	      p(i,j)=pnew
 40	    Continue
 30	  Continue
          write(9,*)iter,dpmax
 20	Continue	    
c
c write data gnuplot 3D format
	Do 50 i=1, maxx
	  Do 60 j=1, maxy
	    Write (8,22) p(i,j)
 60	  Continue
 	Write (8,22)
 50	Continue
 22 	Format(f10.6)
 	Close(8)   
 	Close(9)   
 	Stop 'data saved in laplace.dat'
  	End