IMPLICIT REAL*8 (A-H,O-Z)
      REAL*4 CPUT
      PARAMETER(MJ=1,M=MJ*10,M2=MJ*11-1,N=M*M2,EPS=1.0D-7)
      DIMENSION AR(-M:M,N),AA(N),AB(-M:N),AC(-M:N),F(N), WK(M)
     +,         AB1(-M:N+M),AC1(-M:N+M)
     +,   FQ(N+M)
C
      READ(5,*) DF
      WRITE(6,*) ' % DF= ',DF
      WRITE(6,*) ' % M= ',M,' M2= ',M2,' N= ',N
      FMJ= DFLOAT(MJ)
      DHX=1.0D0/FMJ
      DHY=1.0D0/FMJ
CC
      C0=-1.0D0
      WRITE(6,*) ' % C0= ',C0
      DH=1.0D0/FMJ
      C0H=C0*DH
      C0H2=0.5D0*C0H
C
      IF(DF .EQ. 0.0D0) THEN
        C0H2=0.0D0
      ENDIF
      DO 10 I=1,M-1
        AA(I)= 0.5D0*(DF*4.0D0+4.0D0)
        AB(I)=-0.5D0*((1.0D0-C0H2)+DF*(1.0D0-C0H2))
        AC(I)=-0.5D0*2.0D0
        AB1(I)=-0.5D0*(DF*(1.0D0+C0H2)+(1.0D0+C0H2))
        AC1(I)=0.0D0
   10 CONTINUE     
      AB1(1)=0.0D0
*
      AA(M)= 0.5D0*(DF*(2.0D0-C0H2)+(2.0D0-C0H2))
      AB(M)=0.0D0
      AC(M)=-0.5D0*1.0D0
      AB1(M)=-0.5D0*(DF*(1.0D0+C0H2)+(1.0D0+C0H2))
      AC1(M)=0.0D0
C
      C0H2=0.5D0*C0H
      DO 30 J=1,M2-2
        DO 20 K=1,M-1
          AA(M*J+K)=4.0D0
          AB(M*J+K)=-0.5D0*(2.0D0-C0H)
          AC(M*J+K)=-0.5D0*2.0D0
          AB1(M*J+K)=-0.5D0*(2.0D0+C0H)
          AC1(M*J+K)=-0.5D0*2.0D0
   20   CONTINUE
        AB1(M*J+1)=0.0D0
*
        AA(M*J+M)= 0.5D0*(4.0D0-C0H)
        AB(M*J+M)=0.0D0
        AC(M*J+M)=-0.5D0*1.0D0
        AB1(M*J+M)=-0.5D0*(2.0D0+C0H)
        AC1(M*J+M)=-0.5D0*1.0D0
   30 CONTINUE
C
      IF(DF .EQ. 0.0D0) THEN
        C0H2=0.0D0
      ENDIF
      DO 40 I=N-M+1, N-1
        AA(I)= 0.5D0*(DF*4.0D0+4.0D0)
        AB(I)=-0.5D0*(DF*(1.0D0-C0H2)+(1.0D0-C0H2))
        AC(I)=0.0D0
        AB1(I)=-0.5D0*((1.0D0+C0H2)+DF*(1.0D0+C0H2))
        AC1(I)=-0.5D0*2.0D0
   40 CONTINUE
      AB1(N-M+1)=0.0D0
*

      AA(N)= 0.5D0*(DF*(2.0D0-C0H2)+(2.0D0-C0H2))
      AB(N)=0.0D0
      AC(N)=0.0D0
      AB1(N)=-0.5D0*((1.0D0+C0H2)+DF*(1.0D0+C0H2))
      AC1(N)=-0.5D0*1.0D0
C
      DO 50 I=1,N
        AR(0,I)=AA(I)
        AR(1,I)=AB(I)
        AR(M,I)=AC(I)
*
        AR(-1,I)=AB1(I)
        AR(-M,I)=AC1(I)
   50 CONTINUE
C
*C    Uhen *** 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     DO 65 J= 5*MJ, 6*MJ-1
C       DO 75 I=1,2*MJ
C         FQ(M*J +2*MJ+I)= 0.4D0
C         FQ(M*J +6*MJ+I)=-0.4D0
C  75   CONTINUE
C  65 CONTINUE
CCC
C     DO 85 J= 5*MJ-1, 6*MJ-1
C       DO 95 I=0,2*MJ
C        KP=M*J +2*MJ+I
C        F(KP)=0.25D0*(DHX*DHY)*(FQ(KP)+FQ(KP+M)+FQ(KP+1)+FQ(KP+1+M))
C        KM=M*J +6*MJ+I
C        F(KM)=0.25D0*(DHX*DHY)*(FQ(KM)+FQ(KM+M)+FQ(KM+1)+FQ(KM+1+M))
C  95   CONTINUE
C  85 CONTINUE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DO 60 J= 5*MJ-1, 6*MJ-1
        DO 70 I=0,2*MJ
          F(M*J +2*MJ+I)= 0.2D0*(DHX*DHY)
          F(M*J +6*MJ+I)=-0.2D0*(DHX*DHY)
   70   CONTINUE
   60 CONTINUE
C
      CALL CLOCK0
C
*C    AR*U=F wo Toku; Kotae --> F ***
      CALL BGLU( N, M, AR, EPS, IR, WK )
      CALL BGSLV( N, M, AR, F )
C
      CALL CLOCK(CPUT)
      UMIN=F(1)
      JUMIN=1
      KUMIN=1
      UMAX=F(1)
      JUMAX=1
      KUMAX=1
      DO 130 J=1,M2
        DO 120 K=1,M
          I=M*(J-1)+K
          IF(F(I).LT.UMIN) THEN
            UMIN=F(I)
            JUMIN=J
            KUMIN=K
          ENDIF
          IF(F(I).GT.UMAX) THEN
            UMAX=F(I)
            JUMAX=J
            KUMAX=K
          ENDIF
  120   CONTINUE
  130 CONTINUE
C
      WRITE(6,*) ' % MJ= ',MJ
      WRITE(6,*) ' % CPUT= ',CPUT
      WRITE(6,*) ' % UMIN,(x,y)= ',UMIN,'(',JUMIN,',',KUMIN,')'
      WRITE(6,*) ' % UMAX,(x,y)= ',UMAX,'(',JUMAX,',',KUMAX,')'
C
      WRITE(6,*) ' z=[ '
C
CC    DO 80 K= M,MJ,-MJ
CC      WRITE(6,2000) (F(M*J+K),J=MJ-1,MJ*10-1,MJ),' ; '
CC 80 CONTINUE
*C    MATLAB You ***
      DO 80 K= MJ,M,MJ
        WRITE(6,2000) (F(M*J+K),J=MJ-1,MJ*10-1,MJ),' ; '
   80 CONTINUE
 2000 FORMAT(1H ,10F6.2,A3)
C
      WRITE(6,*) ' ]; '
      WRITE(6,*) ' contour(0:1:9, 0:1:9, z,10);'
C
      STOP   
      END
CCC
      SUBROUTINE BGLU ( N, M, AR, EPS, IR, WK )
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION AR(-M:M,N), WK(M)
C                           *****
C forward elimination for A
      IR = 0
      DO 100 K = 1, N
C
         IF( ABS(AR(0,K)).LE.EPS )  THEN
            IR = IR+1
            RETURN
         END IF
C
C **************************************
         DO 125 J = K+1, MIN(K+M,N)
  125      WK(J-K) = AR(J-K,K)
C          *******
C **************************************
C
         DO 130 I = K+1, MIN(K+M,N)
           AR(K-I,I) = -AR(K-I,I)/AR(0,K)
           T = AR(K-I,I)
           DO 140 J = K+1, MIN(K+M,N)
  140        AR(J-I,I) = AR(J-I,I)+T*WK(J-K)
C                                    *******
  130    CONTINUE
C
  100 CONTINUE
      RETURN
      END
C
      SUBROUTINE BGSLV ( N, M, AR, B )
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION AR(-M:M,N), B(N)
C                           ****
C forward elimination for B
      DO 100 K = 1, N
C
        T = B(K)
        DO 110 I = K+1, MIN(K+M,N)
  110     B(I) = B(I)+AR(K-I,I)*T
  100 CONTINUE
C
C backward substitution for B
      DO 200 K = N, 1, -1
        S = B(K)
        DO 210 J = K+1, MIN(K+M,N)
  210     S = S-AR(J-K,K)*B(J)
        B(K) = S/AR(0,K)
  200 CONTINUE
      RETURN
      END
CCC
      SUBROUTINE CLOCK0
      DIMENSION IA(3)
      COMMON /CLO/ IB(3)
C
        CALL ITIME(IA)
        DO 200 I=1,3
  200     IB(I)=IA(I)
CC        WRITE(6,*) (IB(I),I=1,3)
       RETURN
       END
CCC
      SUBROUTINE CLOCK(CPUT)
      DIMENSION IA(3)
      COMMON /CLO/ IB(3)
C
        CALL ITIME(IA)
        IJI= IA(1)-IB(1)
        IFUN=IA(2)-IB(2)
        IBYO=IA(3)-IB(3)
        CPUT=IJI*3600.0 +IFUN*60.0 +IBYO*1.0
CC        WRITE(6,*) (IB(I),I=1,3),(IA(I),I=1,3),CPUT
C
        CALL ITIME(IA)
        DO 200 I=1,3
  200     IB(I)=IA(I)
      RETURN
      END