PROGRAM MAIN
C  ****************************************************************************
C                                                                             *
C  IMPORTANT NOTES                                                            *
C  ***************                                                            *
C                                                                             *
C  1. IWORK doesn't have to be as big as it is set here unless you are        *
C     running at seriously tight tolerances.                                  *
C                                                                             *
C  2. CPU times are measured using the routine DTIME, which is available      *
C     on SUN and Silicon Graphics machines (although it doesn't always        *
C     seem to do what the documentation says it does). For other machines,    *
C     DTIME will have to be replaced.                                         *
C                                                                             *
C  3. SLEUTH itself does not need the array WORKS to be as big as it is.      *
C     It needs an array of length 1136. The big size of WORKS is due to       *
C     the computation of eigenfunctions later on, in the call to SL4EFN.      *
C     SL4EFN requires 27 DOUBLE PRECISION storage allocations for each        *
C     meshpoint used, plus a further 56 DOUBLE PRECISION storage allocations. *
C     We could have written SL4EFN to require very few storage allocations    *
C     but then the eigenfunction computation process would either have been   * 
C     more expensive or else less flexible in the way that the user can       *
C     specify the evaluation points.                                          *
C                                                                             *
C                                                                             *
C  END IMPORTANT NOTES                                                        *
C  ****************************************************************************
C     .. Parameters ..
      INTEGER IWORK,NCOARS,IWK,ISMAL
      PARAMETER (IWORK=40960,NCOARS=40,IWK=100000,ISMAL=400)
C
      INTEGER IWORKS
      PARAMETER (IWORKS=27*IWORK+56)
      INTEGER YWK
      PARAMETER (YWK=400)
C     ..
C     .. Scalars in Common ..
      INTEGER IPROB
      DOUBLE PRECISION ALFA,BETA,RPARM,M,EPSIL
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION A,B,ELAM,ERR,TOL,X,SUM,DIFF,MAXSUM,MAXDIF
      REAL ELAPSE
      INTEGER I,IFAIL,IMATCH,K,NMESH,NXTRAP,NPLOT
      LOGICAL HOT,SYM
      CHARACTER*1 CYN,AINFO,BINFO
C     ..
C     .. Local Arrays ..
      DOUBLE PRECISION TARRAY(2),WORK(0:IWORK,1:6),WORKS(IWORKS),
     &                 Y(1:4),YSTOR1(4,0:IWORK),YSTOR2(4,0:IWORK),
     &                 Y1(1:4),Y2(1:4)
      INTEGER KNOBS(2)
C     ..
C     .. External Functions ..
      REAL DTIME
      EXTERNAL DTIME
C     ..
C     .. External Subroutines ..
      EXTERNAL EFN,SL4BCS,SL4COF,SLEUTH
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,ATAN,MAX
C     ..
C     .. Common Blocks ..
      COMMON /COMBLK/IPROB
      COMMON /PARM/ALFA,RPARM,M,BETA,EPSIL
C     ..
      WRITE (6,FMT=*) 'Specify problem number:'
      READ (5,FMT=*) IPROB

      IF (IPROB.EQ.1) THEN
          WRITE (6,FMT=*) 'Specify parameter BETA:'
          READ (5,FMT=*) BETA
      END IF

      IF (IPROB.EQ.7) THEN
          WRITE (6,FMT=*) 'Specify parameter ALFA:'
          READ (5,FMT=*) ALFA
      END IF

      IF (IPROB.EQ.8) THEN
          WRITE (6,FMT=*) 'Specify parameter R:'
          READ (5,FMT=*) RPARM
      END IF

      IF (IPROB.EQ.9) THEN
          WRITE (6,FMT=*) 'Specify parameters ALFA and M:'
          READ (5,FMT=*) ALFA,M
      END IF

      IF (IPROB.EQ.20) THEN
          WRITE (6,FMT=*) 'Specify BETA:'
          READ (5,FMT=*) BETA
      END IF

      IF (IPROB.EQ.23) THEN
          WRITE (6,FMT=*) 'Specify EPSILON:'
          READ (5,FMT=*) EPSIL
      END IF

      A = -2.D0*ATAN(1.D0)
      B = -A
      WRITE (6,FMT=*) 'Specify endpoints A and B:'
      READ (5,FMT=*) A,B
      WRITE (6,FMT=*) 'Specify eigenvalue index:'
      READ (5,FMT=*) K
      WRITE (6,FMT=*) 'Is the problem symmetric (T or F):'
      READ (5,FMT=9000) CYN

9000  FORMAT (A)

      SYM = .FALSE.
      IF (CYN.EQ.'T') SYM = .TRUE.
      WRITE (6,FMT=*) 'Specify tolerance:'
      READ (5,FMT=*) TOL

      ELAPSE = DTIME(TARRAY)
      ELAM = 0.D0
      WRITE (6,FMT=*) 'Specify initial lambda:'
      READ (5,FMT=*) ELAM
      ERR = 0.1D0*MAX(1.D0,ABS(ELAM))
      WRITE (6,FMT=*) 'Specify search parameter:'
      READ (5,FMT=*) ERR
      IFAIL = 0

      CALL SLEUTH(A,B,ELAM,ERR,K,SYM,SL4COF,SL4BCS,TOL,NMESH,IMATCH,
     &            NXTRAP,WORK,IWORK,WORKS,IFAIL)

      ELAPSE = DTIME(TARRAY) - ELAPSE

      IF (IFAIL.EQ.0.D0) THEN
          WRITE (6,FMT=*) 'Successful exit from SLEUTH'
          WRITE (6,FMT=*) 'Eigenvalue approximation:',ELAM
          WRITE (6,FMT=*) 'Error estimate:',ERR
          WRITE (6,FMT=*) 'Number of extrapolations:',NXTRAP
          WRITE (6,FMT=*) 'CPU Time:',ELAPSE
          WRITE (6,FMT=*) 'Number of meshpoints:',NMESH
          WRITE (6,FMT=*) 'Matchpoint:',IMATCH

      ELSE
          WRITE (6,FMT=*) 'Failure exit from SLEUTH, IFAIL=',IFAIL
          WRITE (6,FMT=*) 'Eigenvalue approximation:',ELAM
          WRITE (6,FMT=*) 'Error estimate:',ERR
          WRITE (6,FMT=*) 'Number of extrapolations:',NXTRAP
          WRITE (6,FMT=*) 'CPU Time:',ELAPSE
          STOP
      END IF


      HOT = .FALSE.
      WRITE (6,FMT=*) 'Specify number of plotting points'
      READ (5,FMT=*) NPLOT
      WRITE (7,FMT=*) 'Eigenfunction at points of a uniform mesh:'
      DO 10 I = 0,NPLOT
          X = WORK(0,1) + DBLE(I)*(WORK(NMESH,1)-WORK(0,1))/NPLOT  
          X = MAX(X,WORK(0,1))
          X = MIN(X,WORK(NMESH,1))      
          CALL SL4EFN(X,Y1,Y2,HOT,WORK,IWORK,NMESH,IMATCH,ELAM,MULTI,
     &            WORKS,YSTOR1,YSTOR2,SL4BCS,IFAIL)

          IF (IFAIL.NE.0) THEN
              WRITE (6,FMT=*) 'Failure exit from SL4EFN, IFAIL =',IFAIL
              STOP
          END IF
          IF (I.EQ.0) WRITE (6,FMT=*) 'Multiplicity was',MULTI
          WRITE (7,FMT=*) X,Y1(1)
10    CONTINUE

      STOP

      END

C -------------------------------------------------------------------------

      SUBROUTINE SL4COF(X,P,S,Q,W)

C     .. Scalars Arguments ..
      DOUBLE PRECISION X,P,S,Q,W
C     .. Scalars in Common ..
      INTEGER IPROB
      DOUBLE PRECISION ALFA,BETA,RPARM,M,EPSIL
C     .. Common Blocks ..
      COMMON /COMBLK/IPROB
      COMMON /PARM/ALFA,RPARM,M,BETA,EPSIL

      GOTO(10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160,170,
     &     180,190,200,210,220,230,240) IPROB

C Coffey Evans squared
10    W = 1.D0
      P = 1.D0    
      S = 2.D0 * ((BETA*SIN(2.D0*X))**2)-4.D0*BETA*COS(2.D0*X)
      Q = ((BETA**2)*(SIN(2.D0*X)**2)-2.D0*BETA*COS(2.D0*X))**2
     &   -8.D0*(BETA**2)*((COS(2.D0*X)**2)-(SIN(2.D0*X)**2))
     &   -8.D0*BETA*COS(2.D0*X)
      RETURN

C Made-up-problem (oscillatory q(x)) squared
20    P = 1.D0     
      W = 1.D0
      S = 2.D0*COS(X)+4.D0*COS(2.D0*X)+6.D0*COS(3.D0*X)
      Q = (COS(X)+2.D0*COS(2.D0*X)+3.D0*COS(3.D0*X))**2 +COS(X)+
     &    8.D0*COS(2.D0*X)+27.D0*COS(3.D0*X)
      RETURN

C Modified harmonic oscillator squared
30    P = 1.D0
      W = 1.D0    
      S = 2.D0*((X**2)+(X**4))
      Q = ((X**2)+(X**4))**2 - 2.D0 - 12.D0*(X**2)
      RETURN

C Harmonic oscillator squared
40    P = 1.D0
      W = 1.D0
      S = 2.D0*(X**2)
      Q = x**4-2.D0
      RETURN

C Legendre equation squared
50    P = (1-X**2.D0)**2.D0
      W = 1.D0
      Q = 0.D0
      S = 2.D0*(1-X**2.D0)
      RETURN

C Bessel equation squared
60    P = X
      S = 9.D0/X
      Q = 0.D0
      W = X
      RETURN 

C Littlejohn Problem 1
70    W = 1.D0
      P = 1.D0-X*X
      P = P*P
      S = 4.D0*ALFA*(1.D0-X*X)+8.D0
      Q = 0.D0
      RETURN

C Littlejohn Problem 2
80    W = EXP(-X)
      P = X*X*W
      S = ((2.D0*RPARM+2.D0)*X+2.D0)*W
      Q = 0.D0
      RETURN

C Littlejohn Problem 3
90    P = X*X*((1.D0-X)**(2.D0+ALFA))
      S = 2.D0*((1.D0-X)**(1.D0+ALFA))*((ALFA+M+1.D0)*X+1.D0)
      Q = 0.D0
      W = (1.D0-X)**ALFA
      RETURN

C Littlejohn Problem 1 with ALFA = 0
100   W = 1.D0
      P = 1.D0-X*X
      P = P*P
      S = 8.D0
      Q = 0.D0
      RETURN

C Additional Problem 1
110   P = 1.D0
      S = X+X
      Q = X*X
      W = 1.D0
      RETURN

C Additional Problem 2
120   P = 1.D0
      S = X+X
      Q = X*X
      W = 1.D0
      RETURN

C Additional Problem 3
130   P = 1.D0
      S = 0.D0
      Q = X*X
      W = 1.D0
      RETURN

C Additional Problem 4
140   P = 1.D0
      S = 0.D0
      Q = -1.D0
      W = EXP(-3.D0*X)
      RETURN

C Additional Problem 5
150   P = (1.D0+X)**4
      S = -1.25D0*(X+1.D0)**2
      Q = 0.D0
      W = (X+1.D0)**(-1.5D0)
      RETURN

C Additional Problem 6
160   P = 1.D0
      S = 0.2D0
      Q = 0.0064D0
      W = EXP(-3.D0*0.2D0*X)
      RETURN

C Additional Problem 7
170   P = 1.D0
      S = 2.D0*(X**2)
      Q = X**4-2.D0
      W = 1.D0
      RETURN

C Chanane boundary layer problem
180   P = 1.D0
      S = 0.02D0*X**2
      Q=0.0001D0*x**4-0.02D0
      W = 1.D0
      RETURN

C Easy problem
190   P = 1.D0
      S = 0.D0
      Q = 1.D0
      W = 1.D0
      RETURN

C Made up problem

200   P = 1.D0
      S = 2.D0*((X*(1.D0-X))**2 + BETA*X**2*(1.D0-(1.D0-X)**2))
      Q = ((X*(1.D0-X))**2  + BETA*X**2*(1.D0-(1.D0-X)**2))**2
     &  -2.D0*(1.D0-X)**2 + 8.D0*X*(1.D0-X) - 2.D0*X**2 
     &  - 2.D0*BETA*(1.D0-(1.D0-X)**2)
     &  -4.D0*BETA*X*(-2.D0*X+2.D0) + 2.D0*BETA*X**2
      W = 1.D0
      RETURN

210   P = 1.D0
      S = 2.D0*((X*(1.D0-X))**2 + BETA*X**2*(1.D0-(1.D0-X)**2))
      Q = ((X*(1.D0-X))**2  + BETA*X**2*(1.D0-(1.D0-X)**2))**2
     &  -2.D0*(1.D0-X)**2 + 8.D0*X*(1.D0-X) - 2.D0*X**2 
     &  - 2.D0*BETA*(1.D0-(1.D0-X)**2)
     &  -4.D0*BETA*X*(-2.D0*X+2.D0) + 2.D0*BETA*X**2
      W = 1.D0
      RETURN

220   P = 1.D0
      S = 0.D0
      Q = X 
      W = 1.D0
      RETURN

230   P = 1.D0
C      S = 2.D0*X*X/(1.D0+EPSIL*X*X)
C      IF (ABS(X).LT.20.D0) THEN
C       S = 2.D0*X*X
C      ELSE
C       S = 0.D0
C      END IF
      S = -2.D0*X
      Q = X**4-2.D0 
      W = 1.D0
      RETURN

240   P = 1.D0
      S = (1.D0/X)**2
      Q = -20.D0*EXP(-X)+1.D0
      W = 1.D0
      RETURN
      END


      SUBROUTINE SL4BCS(IEND,ISING,X,U,V,ELAM)
C     .. Scalar Arguments ..
      DOUBLE PRECISION ELAM,X
      INTEGER IEND
      LOGICAL ISING
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION U(2,2),V(2,2)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION ALPHA,EPS
C     ..
C     .. Scalars in Common ..
      INTEGER IPROB
      DOUBLE PRECISION ALFA,BETA,RPARM,M,EPSIL
      DOUBLE PRECISION E,E2,CH,CS,SH,SS,NU,OMEGA
C     ..
C     .. Common blocks ..
      COMMON /COMBLK/IPROB
      COMMON /PARM/ALFA,RPARM,M,BETA,EPSIL

      GO TO (10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160,
     &       170,180,190,200,210,220,230,240) IPROB


10        U(1,1) = 0.D0
          U(1,2) = 0.D0
          U(2,1) = 1.D0
          U(2,2) = 0.D0
          V(1,1) = 1.D0
          V(1,2) = 1.D0
          V(2,1) = 0.D0
          V(2,2) = 0.D0
          RETURN 
20        U(1,1) = 0.D0
          U(1,2) = 0.D0 
          U(2,1) = 1.D0
          U(2,2) = 0.D0
          V(1,1) = 0.D0
          V(1,2) = -1.D0
          V(2,1) = 0.D0
          V(2,2) = 0.D0

          RETURN

30        U(1,1) = 0.D0
          U(1,2) = 0.D0
          U(2,1) = 1.D0
          U(2,2) = 0.D0
          V(1,1) = 1.D0
          V(1,2) = 1.D0
          V(2,1) = 0.D0
          V(2,2) = 0.D0
      RETURN

40        U(1,1) = 0.D0
          U(1,2) = 0.D0
          U(2,1) = 0.D0
          U(2,2) = 0.D0
          V(1,1) = 1.D0
          V(1,2) = 0.D0
          V(2,1) = 0.D0
          V(2,2) = 1.D0
      RETURN

C40        U(1,1) = 1.D0
C          U(1,2) = 0.D0
C          U(2,1) = 0.D0
C          U(2,2) = 1.D0
C          V(1,1) = 0.D0
C          V(1,2) = 0.D0
C          V(2,1) = 0.D0
C          V(2,2) = 0.D0
C      RETURN
 
50        U(1,1) = 1.D0
          U(1,2) = 0.D0
          U(2,1) = 0.D0
          U(2,2) = 0.D0
          V(1,1) = 0.D0
          V(1,2) = 0.D0
          V(2,1) = 0.D0
          V(2,2) = 1.D0
      RETURN
60    IF (IEND.EQ.0) THEN
          ALPHA = (1.D0/X)*(2.D0+ELAM*X*X)/(1.D0+(ELAM/4.D0)*X**2.D0)
C          WRITE (10,FMT=*) 'EP = ',X,' AND ALPHA = ',ALPHA
C          ALPHA = 2.D0/X
          U(1,1) = 1.D0
          U(1,2) = 0.D0
          U(2,1) = ALPHA
          U(2,2) = 0.D0
          V(1,1) = 0.D0
          V(1,2) = 1.D0
          V(2,1) = (1.D0/ALPHA)*((8.D0/X)+8.D0*ALPHA-(ALPHA**2.D0))
          V(2,2) = -1.D0/ALPHA
      ELSE
          U(1,1) = 0.D0
          U(1,2) = 0.D0
          U(2,1) = 1.D0
          U(2,2) = 0.D0
          V(1,1) = 0.D0
          V(1,2) = 1.D0
          V(2,1) = -1.D0
          V(2,2) = 0.D0
      END IF
      RETURN

C Lambda-dependent BC for Littlejohn Problem 1
70    U(1,1) = 1.D0
      U(2,1) = ELAM/(8.D0*ALFA*X)
      U(1,2) = 0.D0
      U(2,2) = 0.D0
      V(1,1) = ELAM/(ALFA*X)
      V(2,1) = 0.D0
      V(1,2) = -ELAM/(8.D0*ALFA*X)
      V(2,2) = 1.D0
      RETURN

C Lambda-dependent BC for Littlejohn Problem 2
80    IF (IEND.EQ.0) THEN
C          U(1,1) = 1.D0
C          U(2,1) = -0.5D0*ELAM/RPARM
C          U(1,2) = 0.D0
C          U(2,2) = 0.D0
C          V(1,1) = -ELAM/RPARM
C          V(2,1) = 0.D0
C          V(1,2) = 0.5D0*ELAM/RPARM
C          V(2,2) = 1.D0
         U(1,1) = 1.D0
         U(2,1) = 0.D0
         U(1,2) = -3.D0+3.D0*X*X/10.D0
         U(2,2) = 3.D0*X/10.D0
         V(1,1) = 0.D0
         V(2,1) = 0.D0
         V(1,2) = 3.D0*X*X*EXP(-X)
         V(2,2) = 3.D0*X*X*EXP(-X)

      ELSE
          U(1,1) = 1.D0
          U(1,2) = 0.D0
          U(2,1) = 0.D0
          U(2,2) = 1.D0
          V(1,1) = 0.D0
          V(2,1) = 0.D0
          V(1,2) = 0.D0
          V(2,2) = 0.D0
      END IF
      RETURN

C Lambda-dependent BC for Littlejohn Problem 3
90    IF (IEND.EQ.0) THEN
         U(1,1) = 1.D0
         U(2,1) = -0.5D0*ELAM/M
         U(1,2) = 0.D0
         U(2,2) = 0.D0

         V(1,1) = -ELAM/M
         V(2,1) = 0.D0
         V(1,2) = 0.5D0*ELAM/M
         V(2,2) = 1.D0

      ELSE
         U(1,1) = 1.D0
         U(2,1) = 0.D0
         U(1,2) = 0.D0
         U(2,2) = 1.D0

         V(1,1) = 0.D0
         V(2,1) = 0.D0
         V(1,2) = 0.D0
         V(2,2) = 0.D0         
      END IF

      RETURN

C Lambda-independent BC for Littlejohn Problem 1 with ALFA = 0.
100   U(1,1) = 0.D0
      U(2,1) = 1.D0
      U(1,2) = 0.D0
      U(2,2) = 0.D0
      V(1,1) = 0.D0
      V(2,1) = 0.D0
      V(1,2) = 1.D0
      V(2,2) = 0.D0

      RETURN

C Additional problem 1
110   U(1,1) = 1.D0
      U(1,2) = 0.D0
      U(2,1) = 0.D0
      U(2,2) = 1.D0
      V(1,1) = 0.D0
      V(2,1) = 0.D0
      V(1,2) = 0.D0
      V(2,2) = 0.D0
      RETURN

C Additional problem 2
120   U(1,1) = 0.D0
      U(1,2) = 0.D0
      U(2,1) = 0.D0
      U(2,2) = 0.D0
      V(1,1) = 1.D0
      V(2,1) = 0.D0
      V(1,2) = 0.D0
      V(2,2) = 1.D0
      RETURN

C Additional problem 3
130   U(1,1) = 0.D0
      U(1,2) = 0.D0
      U(2,1) = 0.D0
      U(2,2) = 0.D0
      V(1,1) = 1.D0
      V(2,1) = 0.D0
      V(1,2) = 0.D0
      V(2,2) = 1.D0
      RETURN

C Additional problem 4
140   NU = 1.D0
      OMEGA = 1.D0
      E = EXP(-NU*X)
      E2 = E*E
      CH = (1.D0+E2)/2.D0
      SH = (1.D0-E2)/2.D0
      CS = COS(OMEGA*X)*E
      SS = SIN(OMEGA*X)*E

      U(1,1) = CS + CH
      U(1,2) = SS + SH
      U(2,1) = -SS + SH
      U(2,2) = CS + CH

      V(1,1) = -SS - SH
      V(1,2) = CS - SH
      V(2,1) = -CS + CH
      V(2,2) = -SS + SH

C Additional Problem 5
150   IF (IEND.EQ.0) THEN
C Neumann conditions
      U(1,1) = 1.D0
      U(2,1) = 0.D0
      U(1,2) = 0.D0
      U(2,2) = 1.D0

      V(1,1) = 0.D0
      V(2,1) = 0.D0
      V(1,2) = 0.D0
      V(2,2) = 0.D0

      ELSE       

C Friedrichs conditions
      U(1,1) = 1.D0/(1.D0+X)
      U(2,1) = -1.D0/(1.D0+X)**2
      V(1,1) = -0.75D0
      V(2,1) = 2.D0*(1.D0+X)
     
C      U(1,2) = 0.D0
C      U(2,2) = 0.D0
C      V(1,2) = 1.D0/(1.D0+X)**2
C      V(2,2) = 1.D0/(1.D0+X)
      U(1,2) = 1.D0/(X+1.D0)**1.5D0
      U(2,2) = -1.5D0*U(1,2)/(1.D0+X)
      V(1,2) = 0.D0
      V(2,2) = 15.D0*SQRT(X+1.D0)/4.D0

      END IF

      RETURN

C Additional Problem 6
160   IF (IEND.EQ.0) THEN
C Neumann conditions
      U(1,1) = 1.D0
      U(2,1) = 0.D0
      U(1,2) = 0.D0
      U(2,2) = 1.D0

      V(1,1) = 0.D0
      V(2,1) = 0.D0
      V(1,2) = 0.D0
      V(2,2) = 0.D0

      ELSE       
C Friedrichs conditions
      U(1,1) = 1.D0
      U(2,1) = 0.2D0
      V(1,1) = 0.032D0
      V(2,1) = 0.04D0

C      U(1,2) = 0.0D0
C      U(2,2) = 0.0D0
C      V(1,2) = -0.2D0
C      V(2,2) = 1.D0

      U(1,2) = 1.D0
      U(2,2) = -0.4D0
      V(1,2) = -0.016D0
      V(2,2) = 0.16D0
     
      END IF

      RETURN

C Additional Problem 7
170   U(1,1) = 0.D0
      U(1,2) = 0.D0
      U(2,1) = 0.D0
      U(2,2) = 0.D0
      V(1,1) = 1.D0
      V(1,2) = 0.D0
      V(2,1) = 0.D0
      V(2,2) = 1.D0

      RETURN

C Chanane boundary layer problem
180   U(1,1) = 0.D0
      U(1,2) = 0.D0
      U(2,1) = 0.D0
      U(2,2) = 1.D0
      V(1,1) = 1.D0
      V(1,2) = 0.D0
      V(2,1) = 0.D0
      V(2,2) = 0.D0

      RETURN

C Easy Problem
190   U(1,1) = 0.D0
      U(1,2) = 0.D0
      U(2,1) = 0.D0
      U(2,2) = 0.D0
      V(1,1) = 1.D0
      V(1,2) = 0.D0
      V(2,1) = 0.D0
      V(2,2) = 1.D0

      RETURN

C Next problem has Neumann conditions at both ends:
200   U(1,1) = 1.D0
      U(2,1) = 0.D0
      U(1,2) = 0.D0   
      U(2,2) = 1.D0
      V(1,1) = 0.D0
      V(2,1) = 0.D0
      V(1,2) = 0.D0   
      V(2,2) = 0.D0
      RETURN
C Next problem has Neumann conditions at right hand end and
C Dirichlet condition at left hand end:
210   IF (IEND.EQ.1) THEN
          U(1,1) = 1.D0
          U(2,1) = 0.D0
          U(1,2) = 0.D0   
          U(2,2) = 1.D0
          V(1,1) = 0.D0
          V(2,1) = 0.D0
          V(1,2) = 0.D0   
          V(2,2) = 0.D0
      ELSE
          U(1,1) = 0.D0
          U(2,1) = 0.D0
          U(1,2) = 0.D0   
          U(2,2) = 0.D0
          V(1,1) = 1.D0
          V(2,1) = 0.D0
          V(1,2) = 0.D0   
          V(2,2) = 1.D0
      END IF
      RETURN

220   U(1,1) = 0.D0
      U(2,1) = 0.D0
      U(1,2) = 0.D0
      U(2,2) = 0.D0
      V(1,1) = 1.D0
      V(1,2) = 0.D0
      V(2,1) = 0.D0
      V(2,2) = 1.D0
      RETURN

230   U(1,1) = 0.D0
      U(2,1) = 0.D0
      U(1,2) = 0.D0
      U(2,2) = 0.D0
      V(1,1) = 1.D0
      V(1,2) = 0.D0
      V(2,1) = 0.D0
      V(2,2) = 1.D0
      RETURN

240   U(1,1) = 0.D0
      U(2,1) = 0.D0
      U(1,2) = 0.D0
      U(2,2) = 0.D0
      V(1,1) = 1.D0
      V(1,2) = 0.D0
      V(2,1) = 0.D0
      V(2,2) = 1.D0
      RETURN

      END