CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      PROGRAM BOUND_UNI 
C     Bound states in a one-dimensional potential
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
5     CONTINUE            !main loop; execute once for each set of param
         CALL PARAM       !get input parameters
         CALL ARCHON      !solve time-independent Schroedinger equation
      GOTO 5
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE ARCHON
C solves the time-independent Schroedinger equation for 1-dimensional
C potential; allows for more than one eigenvalue search
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PAR.E3'
C Local variables:
      REAL ENERGY(MAXLEV)            !array of eigenvalues
      INTEGER NFOUND                 !number of levels found
      REAL PSI(MAXLAT)               !wavefunction
      INTEGER ITER                   !continue finding levels?
      INTEGER NODES                  !number of nodes
      REAL E                         !trial eigenenergy 
      INTEGER LTP,RTP                !classical turning points (lattice)
      REAL DE,DESAVE                 !step in energy search
      INTEGER NSTP                   !number of steps taken in search
C Global variables:
      INCLUDE 'PAR2.E3'
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      OPEN(OUNIT,FILE='EX31.OUT',STATUS='UNKNOWN')
      OPEN(MUNIT,FILE='EX32.OUT',STATUS='UNKNOWN')
      OPEN(NUNIT,FILE='EX33.OUT',STATUS='UNKNOWN')
C
      ITER=1                          !initialize searching variables
      NFOUND=0
      E=EBOTTM                        !reasonable starting values
      DE=DEBOT
C                                                    
5     IF (ITER .EQ . 1) THEN          !loop over energy levels
C                                      
C          PRINT*,' Enter starting values of E, EMIN, EMAX'
C          PRINT*,' E = approximate eigenvalue'
C          PRINT*,' EMIN = lower limit of searching region' 
C          PRINT*,' EMAX = upper limit of searching region' 
C          READ(*,*)E,EMIN,EMAX
          DESAVE=DE                   !save for next level                                                            
          CALL SEARCH(E,DE,LTP,RTP,NODES,NSTP,PSI) !search for energy
C
C         prompt for new E and DE if too many steps taken in search
          IF (NSTP .GE. MAXSTP) THEN
              PRINT*, ' '
              PRINT*, ' Eigenstate not found - too many steps'
              PRINT*, ' Change value of E and DE, or increase MAXSTP'
          ELSE
C         or output eigenvalue information
          CALL EOUT (E,LTP,RTP,NODES)
C           save energy and turning points
            IF (NFOUND .EQ. MAXLEV) THEN
               NFOUND=0
               PRINT*, ' Storage capacity for levels exceeded'
               PRINT*, ' Old information will be overwritten'
            END IF                        
            NFOUND=NFOUND+1
            ENERGY(NFOUND)=E
C
C           check if another level is to be found
            PRINT*, ' '
            PRINT*, ' '
            PRINT*,'Do you want to find another level? (1=YES)'
            READ(*,*)IOPT
            IF(IOPT.NE.1) THEN
               CALL OUT(PSI)
               PRINT*,
     +       ' ***   Outputs in EX31.OUT, EX32.OUT and EX33.OUT  ***'
               STOP
            ENDIF
C
          END IF
C
          E=E+DESAVE          !reasonable starting values for next level
          DE=DESAVE
C
      GOTO 5
      END IF        
      CALL OUT(PSI)
C
      RETURN
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE SEARCH(E,DE,LTP,RTP,NODES,NSTP,PSI)
C search for one eigenvalue beginning with energy E and step DE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PAR.E3'
C Passed variables:
      INTEGER NODES          !number of nodes (output)
      REAL E                 !trial eigenenergy, true eigenenergy (I/O)
      INTEGER LTP,RTP        !classical turning points (lattice)(output)
      REAL DE                !step in energy search (input)
      INTEGER NSTP           !number of steps taken in search (output)
      INTEGER NCOUNT         !output heading in FOUT
      REAL PSI(MAXLAT)       !wavefunction (output)
C Local variables:
      LOGICAL SECANT         !are we doing secant search yet?
      REAL F,FOLD,FSTART     !values of the mismatch
      REAL EOLD              !last value of the energy
      INTEGER NCROSS         !classically forbidden regions
      LOGICAL FIRST          !signal first step in search
C Global variables:
      INCLUDE 'PAR2.E3'
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      FIRST=.TRUE.           !first call to NUMERV for this level
      SECANT=.FALSE.         !start off with simple search
      F=TOLF                 !dummy value to get into loop
      NSTP=0
      NCOUNT=0
C
C    search for eigenenergy until F is small enough or too many
C    steps are taken
10    IF ((ABS(F) .GE. TOLF) .AND. (NSTP .LT. MAXSTP)) THEN
          NSTP=NSTP+1
C         integ Schroedinger equation; find discontinuity in derivative
          CALL NUMERV (F,NODES,E,LTP,RTP,NCROSS,FIRST,PSI)
          IF (FIRST .EQV. .TRUE.) THEN
              FSTART=F
              FIRST=.FALSE.
          END IF
          IF (F*FSTART  .LE. 0.) SECANT=.TRUE.    
          IF ((SECANT.EQV..TRUE.).AND.(ABS(F-FOLD).GT.0.)) THEN
              DE=-F*(E-EOLD)/(F-FOLD)
          END IF
C
          CALL FOUT(NCOUNT,E,DE,F,NODES,NCROSS)
C
          EOLD=E             !update values for next step
          FOLD=F                                                     
          E=E+DE
          NCOUNT=NCOUNT+1
C                                
C         keep energy greater than VMIN=-1
          IF (E .LT. -1.) THEN
              PRINT*, ' Energy is less than -1.'
          END IF
C         end search if DE is smaller than machine precision         
          IF ((ABS(DE) .LT. DEMIN) .AND. (ABS(F) .GT. TOLF)) THEN
             WRITE (*,20)
20           FORMAT('0',21X,'Delta E is too small, search must stop')
             F=TOLF/2             !dummy value to end search
          END IF
C
      GOTO 10  
      END IF 
C
      RETURN
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE NUMERV (F,NODES,E,LTP,RTP,NCROSS,FIRST,PSI)
C subroutine to integrate time independent Schroedinger equation
C with energy=E
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PAR.E3'
C Passed variables:
      REAL F                     !size of mismatch (output)
      INTEGER NODES              !number of nodes (output)
      REAL E                     !trial eigenvalue (input)
      INTEGER LTP,RTP            !classical turning points (output)
      INTEGER NCROSS             !classically forbidden regions (output)
      LOGICAL FIRST              !signal first step in search (output)
      REAL PSI(MAXLAT)           !wavefunction (output)
C Local variables:
      INTEGER IMATCH             !matching lattice point
      REAL C                     !useful constant
      INTEGER IX,KX              !X indices
      REAL KI,KIM1,KIP1          !terms in Numerov algorithm
      REAL NORM                  !norm of wavefunction
      REAL PMMTCH,PPMTCH         !PSI at match point
      LOGICAL FLIP               !do we need to flip sign of PSI
      Real PSIMAX                !max value of wavefunction
C Global variables:
      INCLUDE 'PAR2.E3'
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      IF (FIRST .EQV. .TRUE.) IMATCH=0   !signals search for IMATCH
      C=(DX*DX/12)*GAMMA*GAMMA           !evaluate constant
C                                                                        
C     find IMATCH once per energy level by looking for entry
C     into classically forbidden regions
      IX=1
      KIM1=C*(E-V(IX))
1     IF (IMATCH .EQ. 0) THEN
         IX=IX+1
         KI=C*(E-V(IX))    !neg value of K indicates class frbdn region
         IF ((KI*KIM1 .LT. 0) .AND. (KIM1 .GT. 0)) IMATCH=IX
C        if other procedure fails to find IMATCH, set IMATCH=NPTS-10
         IF (IX .EQ. NPTS-10) IMATCH=IX
         KIM1=KI
      GOTO 1
      END IF
C
      PSI(1)=0                     !left hand side bound. cond.
      PSI(2)=9.999999E-10          
      KIM1=C*(E-V(1))              !initial K*K values
      KI=C*(E-V(2))
C
C     Numerov algorithm; S=0; integrate until enter class. frbdn. region
      DO 10 IX=2,IMATCH
         KIP1=C*(E-V(IX+1))
         PSI(IX+1)=(PSI(IX)*(2.-10.*KI)-PSI(IX-1)*(1+KIM1))/(1+KIP1)
         KIM1=KI                  !roll values of K*K
         KI=KIP1
C
C        if PSI grows too large rescale all previous points
         IF (ABS(PSI(IX+1)) .GT. (1.0E+10)) THEN
             DO 20 KX=1,IX+1
                PSI(KX)=PSI(KX)*9.999999E-06
20           CONTINUE
         END IF
10    CONTINUE
      PMMTCH=PSI(IMATCH)           !save value for normalization
C
      PSI(NPTS)=0                  !rhs boundary conditions
      PSI(NPTS-1)=9.999999E-10
      KIP1=C*(E-V(NPTS))           !initial K*K values
      KI=C*(E-V(NPTS-1))
C
C     Numerov algorithm, S=0;integrate from rhs to IMATCH
      DO 30 IX=NPTS-1,IMATCH+1,-1   
           KIM1=C*(E-V(IX-1))
           PSI(IX-1)=(PSI(IX)*(2.-10.*KI)-PSI(IX+1)*(1+KIP1))/(1+KIM1)
           KIP1=KI                 !roll values of K*K
           KI=KIM1
C 
C          if PSI grows too large rescale all previous points
           IF (ABS(PSI(IX-1)) .GT. (1.0E+10)) THEN 
              DO 40 KX=NPTS-1,IX-1,-1
                   PSI(KX)=PSI(KX)*9.999999E-06
40            CONTINUE
           END IF
30    CONTINUE                                    
C
      KIM1=C*(E-V(IMATCH-1))       !finds values needed for log deriv
      PPMTCH=(PSI(IMATCH)*(2-10*KI)-PSI(IMATCH+1)*(1+KIP1))/(1+KIM1)
C
      NORM=PMMTCH/PSI(IMATCH)      !norm PSI right to PSI left
      DO 5 IX=IMATCH,NPTS                                     
           PSI(IX)=PSI(IX)*NORM
5     CONTINUE
      PPMTCH=PPMTCH*NORM
C                                                       
C     find nodes, turning points, entry into classically forb regions
C     and determine if overall sign must be flipped
      CALL DETAIL(IMATCH,NODES,E,LTP,RTP,NCROSS,PMMTCH,FLIP,FIRST,PSI
     +    ,nlft,nrt      )
C
C     find maximum PSI value and flip sign if necessary
      PSIMAX=ABS(PSI(1))           
      DO 7 IX=1,NPTS
           IF (ABS(PSI(IX)) .GT. PSIMAX) PSIMAX=ABS(PSI(IX))
           IF (FLIP .EQV. .TRUE.) PSI(IX)=-PSI(IX)
7     CONTINUE
      IF (FLIP) PPMTCH=-PPMTCH
C
      F=(PSI(IMATCH-1)-PPMTCH)/PSIMAX   !evaluate matching condition
C
      RETURN 
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE 
     +    DETAIL(IMATCH,NODES,E,LTP,RTP,NCROSS,PMMTCH,FLIP,FIRST,PSI
     +     ,nlft,nrt)
C subroutine to calculates nodes, turning points, entry
C into classically forbidden regions, and whether or not the
C overal sign of PSI must be flipped to make F continuous function of E
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PAR.E3'
C Passed variables:
      INTEGER IMATCH            !matching lattice point (input)
      INTEGER NODES             !number of nodes (output)
      REAL E                    !trial eigenvalue (input)
      INTEGER LTP,RTP           !classical turning points (output)
      INTEGER NCROSS            !classically forbidden regions (output)
      REAL PMMTCH               !PSI at match point (input)
      LOGICAL FIRST             !signal first step in search(input)
      LOGICAL FLIP              !do we need to flip sign of PSI (output)
      REAL PSI(MAXLAT)          !wavefunction (input)
C Local variables:
      LOGICAL LSAME             !does PSI left retain its sign?
      INTEGER IX                !X index
      REAL K,KLAST              !wavenumbers 
      INTEGER NLOLD,NLFT        !left nodes
      INTEGER NROLD,NRT         !right nodes
      REAL PSIR,PSIL            !sign of PSI which must be constant
C Global variables:
      INCLUDE 'PAR2.E3'
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  initialize variables
      PSIL=PSI(2)
      PSIR=PSI(NPTS-1)
C  find number of nodes on each side
      NLOLD=NLFT                !save old values first
      NROLD=NRT
      NLFT=0                    !zero sums                 
      NRT=0
      DO 8 IX=2,IMATCH-1        !left side
          IF (PSI(IX)*PSI(IX-1) .LT. 0) NLFT=NLFT+1
8     CONTINUE
      IF (PSI(IMATCH-1)*PMMTCH .LT. 0) NLFT=NLFT+1
      DO 10 IX=IMATCH+1,NPTS        !right side
          IF (PSI(IX)*PSI(IX-1) .LT. 0) NRT=NRT+1             
10    CONTINUE
C
C     check for change in node number
C      if one side of the wavefunction has gained a node,
C      that side must maintain the same sign between steps
C      in order that F be a contintuous function of Energy.
C      if no new node has appeared, keep same sign same
C      as last time
      IF ((NLOLD .NE. NLFT) .AND. (NRT .EQ. NROLD)) THEN
         IF (LSAME .EQV. .FALSE.) LSAME=.TRUE.
      ELSE IF ((NROLD .NE. NRT) .AND. (NLOLD .EQ. NLFT)) THEN
         IF (LSAME .EQV. .TRUE.) LSAME=.FALSE.
      END IF
C
      IF (FIRST .EQV. .TRUE.) LSAME=.TRUE.    !initialize variables
      IF (FIRST .EQV. .TRUE.) PSIL=1.
C                                                           
C     now, finally, determine if the sign needs to be flipped
      FLIP=.FALSE.
      IF ((LSAME) .AND. (PSIL*PSI(2) .LT. 0.)) FLIP=.TRUE.
      IF (( .NOT. LSAME ) .AND. (PSIR*PSI(NPTS-1) .LT. 0.))
     +    FLIP=.TRUE.  
C
C     save values for next step
      PSIL=PSI(2)
      PSIR=PSI(NPTS-1)
      IF (FLIP) PSIL=-PSIL
      IF (FLIP) PSIR=-PSIR
C
      NODES=NLFT+NRT             !total number of nodes
C 
C     find leftmost turning point
      LTP=0
      IX=0
110   IF (LTP .EQ. 0) THEN
         IX=IX+1
         IF ((E-V(IX)) .GE. 0) LTP=IX
      GOTO 110
      END IF
C
C     find rightmost turning point
      RTP=0
      IX=NPTS+1
120   IF (RTP .EQ. 0) THEN
         IX=IX-1
         IF ((E-V(IX)) .GE. 0) RTP=IX
      GOTO 120
      END IF
C
C     do we integrate into classically forbidden regions?
      KLAST=(E-V(LTP))
      NCROSS=0
      DO 90 IX=LTP+1,RTP-1
         K=(E-V(IX))
         IF ((K*KLAST .LE. 0) .AND. (K .LT. 0)) NCROSS=NCROSS+1
         KLAST=K
90    CONTINUE
C
      RETURN
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE POTNTL
C sets up array for current potential;
C all potentials have a minimum value of -1 and max of 1;
C if you change the potential, also change XMIN and XMAX in INIT,
C and estimate DEBOT and EBOTTM (one fifth of the expected spacing
C at the bottom of the well, and two DEBOT's below the expected ground
C state energy).
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PAR.E3'
C Local variables:
      INTEGER ILEFT,IWID,IRIGHT     !lattice values for width, left
      REAL VSCALE                   !scaling factor for renorm of parab
      INTEGER IX                    !labels lattice points
      REAL VMIN                     !minimum value of potential 
      REAL CURVE                    !curvature for smoothing bump
C Global variables:
      INCLUDE 'PAR2.E3'
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     define limits on X and DX
      IF (POT .EQ. SQUARE) THEN
         VXMIN=XMIN1
         VXMAX=XMAX1
      ELSE IF (POT .EQ. PARAB) THEN                  
         VXMIN=XMIN2
         VXMAX=XMAX2
      ELSE IF (POT .EQ. LENRD) THEN
         VXMIN=XMIN3
         VXMAX=XMAX3
      END IF           
C
      DX=(VXMAX-VXMIN)/(NPTS-1)
C
C     setup X (space coordinate) array
      DO 5 IX=1,NPTS
         X(IX)=VXMIN+(IX-1)*DX
5     CONTINUE             
C                                
C     setup V (potential) array
      IF (POT .EQ. SQUARE) THEN             !square well
         DO 10 IX=1,NPTS
            V(IX)=-1.
10       CONTINUE
         VMAX=1.
         DEBOT=.2*(3.1415926/(4*GAMMA))**2!energy spacing at well bottom
         IF (BUMP .EQ. 1) THEN            !bump in well bottom
            ILEFT=NINT((LEFT-XMIN1)/DX)+1 !locate bump edges
            IRIGHT=NPTS+NINT((RIGHT-XMAX1)/DX)
            RIGHT=X(IRIGHT)
            LEFT=X(ILEFT)
            CURVE=-.025/((RIGHT-LEFT)/2.)**2
C           top of bump is curved to avoid discontinuities
            DO 20 IX=ILEFT,IRIGHT
               V(IX)=-1.+HEIGHT*(1.+CURVE*(X(IX)-(RIGHT+LEFT)/2.)**2)
            print*,ix,v(iX)
20          CONTINUE
         END IF
C
      ELSE IF (POT .EQ. PARAB) THEN    !parabolic well             
         DO 30 IX=1,NPTS
            V(IX)=-(1.-.5*X(IX)*X(IX))
30       CONTINUE
         VMAX=1. 
         DEBOT=.2/(SQRT(2.)*GAMMA)     !energy spacing at well bottom
         IF (FLAT .EQ. 1) THEN
C           the parabolic potential is flattened symmetrically about x=0
            IWID=INT(PWID/DX)
            IF (MOD(IWID,2) .NE. MOD(NPTS,2)) IWID=IWID-1  
            ILEFT=(NPTS-IWID)/2+1
            VMIN=V(ILEFT)
            DO 40 IX=ILEFT,ILEFT+IWID-1
               V(IX)=VMIN
40          CONTINUE
C
            VSCALE=2./(V(1)-VMIN)    !normalize so that min value is -1
            DO 50 IX=1,NPTS
               V(IX)=VSCALE*(V(IX)-VMIN)-1.
50          CONTINUE
         END IF
C 
      ELSE IF (POT .EQ. LENRD) THEN    !Lennard-Jones potential
         DO 60 IX=1,NPTS  
            V(IX)=4*(X(IX)**(-12.)-X(IX)**(-6))
60       CONTINUE
         VMAX=V(1)
         DEBOT=1.0692/GAMMA !energy spacing at well bottom (see Ex. 1.8)
C
      END IF           
C 
C     educated guesses for ground state energy and energy spacing;
C     with these values, it will take three steps for F to change sign
      EBOTTM=-1.+2.5*DEBOT
C                           
      RETURN                                                 
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE PARAM
C initializes constants, displays header screen,
C initializes arrays for input parameters
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PAR.E3'
      INCLUDE 'PAR2.E3'
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C setup constant values for the potentials
c                             !all potentials have VMIN=-1, VMAX=1
      XMIN1=-2.
      XMAX1=2.
      XMIN2=-2.
      XMAX2=2.
      XMIN3=(2.*(SQRT(2.)-1.))**(.1666666)      
      XMAX3=1.9 
C
C     reasonable limits on energy and energy increment
      EMIN=-.99999
      EMAX=100.
      DEMIN=5.E-06
      DEMAX=1.
C
C     display header screen     
      PRINT*,'EXAMPLE 3'
      PRINT*,'Bound states in a one-dimensional potential'
C
C  potential parameters
      PRINT*,'POTENTIAL FUNCTION MENU'
      PRINT*,' '
      PRINT*,' Enter number (1) [2] {3} for:'
      PRINT*,' (Square-well)'
      PRINT*,' [Parabolic-well]'
      PRINT*,' {Lennard-Jones potential}'
      READ(*,*)POT
C insert bump in the square-well 
      IF(POT.EQ.1) THEN
        PRINT*,' Do you want a bump in the square well?'
        PRINT*,' 0 for no, 1 for yes'
        READ(*,*)BUMP
        IF(BUMP.EQ.1) THEN
          PRINT*,' Enter the left hand edge of bump (scaled units)'
          READ(*,*)LEFT
          PRINT*,' Enter the right hand edge of bump (scaled units)'
          READ(*,*)RIGHT
          PRINT*,' Enter the height of the bump (scaled units)'
          READ(*,*)HEIGHT
        ENDIF
      ENDIF
C  parabolic flattened area
      IF(POT.EQ.2) THEN
        PRINT*,' Do you want a flattened bottom in parabolic-well'
        PRINT*,' 0 for no, 1 for yes'
        READ(*,*)FLAT
        IF(FLAT.EQ.1) THEN
          PRINT*,' Enter length of the flattened area (scaled units)'
          READ(*,*)PWID
        ENDIF
      ENDIF
      CALL POTNTL  !setup potential array
C                
      RETURN
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       SUBROUTINE OUT(PSI)
C outputs results to the specified unit
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
       INCLUDE 'PAR.E3'
C Input variables:
      REAL PSI(MAXLAT)             !wavefunction
C Local variables
      INTEGER IX                   !index of lattice
      REAL NORM,MAXPSI             !norm and max value for PSI
C Global variables:
      INCLUDE 'PAR2.E3'
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C      !message for the impatient
      WRITE (*,100)
C
      NORM=0.
      DO 10 IX=1,NPTS               !normalize wavefunction 
        NORM=NORM+PSI(IX)**2 
10    CONTINUE
      NORM=SQRT(NORM*DX)
      MAXPSI=0.
      DO 30 IX=1,NPTS
         PSI(IX)=PSI(IX)/NORM
         IF (ABS(PSI(IX)) .GT. MAXPSI) MAXPSI = ABS(PSI(IX))
30    CONTINUE
C
       WRITE (OUNIT,2)
       WRITE (OUNIT,4)
       WRITE (OUNIT,6) GAMMA
       WRITE (OUNIT,8) TOLF
       WRITE (OUNIT,12) NPTS
C
       IF (POT .EQ. SQUARE) THEN
          WRITE (OUNIT,14)
          IF (BUMP .EQ. 1) WRITE (OUNIT,16) LEFT,RIGHT,HEIGHT
       ELSE IF (POT .EQ. PARAB) THEN
          WRITE (OUNIT,18)
          IF (FLAT .EQ. 1) WRITE (OUNIT,20) -PWID,PWID
       ELSE IF (POT .EQ. LENRD) THEN
          WRITE (OUNIT,22)
       END IF
C
       WRITE (OUNIT,24)VXMIN,VXMAX,DX
       WRITE (OUNIT,26)
       WRITE (OUNIT,2)
C
2      FORMAT (' ')
4      FORMAT (' Output from Example 3: Time Independent Schroedinger',
     +         ' Equation')
6      FORMAT (' Gamma =',F9.3)
8      FORMAT (' Matching tolerance = ', 1PE10.3)
12     FORMAT (' Number of lattice points =', I5)
14     FORMAT (' Square Well potential')
16     FORMAT (' with bump from ',F6.3,' to ',F6.3,' of height =',F8.3)
18     FORMAT (' Parabolic potential')
20     FORMAT (' with flat bottom from ',F6.3,' to ',F6.3)
22     FORMAT (' Lennard-Jones potential')
24     FORMAT (' Xmin=',F6.3,' Xmax=',F6.3,' Dx=',1PE10.3)
26     FORMAT (' Energies and lengths are in scaled units')
C
      WRITE (OUNIT,79)
      WRITE (OUNIT,80)
      WRITE (OUNIT,85)
      WRITE (OUNIT,70) (X(IX),V(IX),PSI(IX),IX=1,NPTS)
C
70    FORMAT (3(11X,1PE12.5))
79    FORMAT(/,' Wavefunction of the last found eigenstate',/)
80    FORMAT 
     +     ('0',13X,'Position',14X,'Potential',13X,'Wavefunction')
85    FORMAT 
     +     (13X,'---------',14X,'---------',13X,'------------')
100   FORMAT(/,' Patience, please; output going to a file.')
           
C
C
      RETURN
      END
                                                                        
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE EOUT(E,LTP,RTP,NODES)
C output information about eigenstate
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PAR.E3'
C Passed variables:
c      INTEGER MUNIT                 !unit to which we are writing
      INTEGER NODES                 !number of nodes
      REAL E                        !trial eigenenergy 
      INTEGER LTP,RTP               !classical turning points (lattice)
      REAL XLTP,XRTP                !classical turning points
C Global variables:
      INCLUDE 'PAR2.E3'
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      XLTP=VXMIN+(LTP-1)*DX         !turning points
      XRTP=VXMIN+(RTP-1)*DX
      WRITE (MUNIT,500) E, NODES
      WRITE (MUNIT,510) XLTP,XRTP
C
500   FORMAT ('0',20X,'Eigenvalue = ',F10.7,' with ',I2,' nodes.')
510   FORMAT (16X,'and classical turning points ',F7.4,' and ',F7.4)
C
      RETURN
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE FOUT(NCOUNT,E,DE,F,NODES,NCROSS)
C writes results for one level to the requested unit
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PAR.E3'
C Input variables:
C      INTEGER NUNIT                !unit to which we are writing
      REAL F                       !values of the mismatch
      INTEGER NODES                !number of nodes
      REAL E                       !trial eigenenergy 
      REAL DE                      !step in energy search
      INTEGER NCROSS               !classically forbidden regions
      INTEGER NCOUNT               !output heading in FOUT
      CHARACTER*3 FORBDN           !have we integrated into forb regions
C Global variables:
      INCLUDE 'PAR2.E3'
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C
      IF(NCOUNT.EQ.0) THEN
        WRITE (NUNIT,31)
        WRITE (NUNIT,35)
31      FORMAT 
     +     (10X,'Energy',14X,'De',16X,'F',10x,'Nodes',5x,'Frbddn')
35      FORMAT 
     +     (10X,'------',14X,'--',16X,'-',10X,'-----',5X,'------')
      ENDIF
C
      IF (NCROSS .GT. 0) THEN
         FORBDN='Yes'
      ELSE
         FORBDN='No '
      END IF
      WRITE (NUNIT,20) E,DE,F,NODES,FORBDN
C                 
20    FORMAT (8X,F10.7,7X,1PE10.3,7X,1PE10.3,8X,I2,8X,A3)
C
      RETURN
      END