program eigen_parabola
      implicit none
      integer pgopen, i, odd_even, n_steps, i_again, i_colour
      real x, phi, dphi_dx, V0, a, eigenvalue, phi_ini, dphi_dx_ini
      real delta_x, Vx(6), Vy(6), V, phi_new, half_a

      write(*,*) 'Type in depth V0 and width a of the potential well:'
      read (*,*) V0, a
      write(*,*) 'Type in small spacing delta_x:'
      read(*,*) delta_x

      if (pgopen('/xwin').le.0) stop
      call pgpap (5.0,0.75)
      call pgenv (-20.0,20.0,-50.0,20.0,0,0)

c     Vx(1) = -20.0
c     Vy(1) = 0.0
c     Vx(2) = -1*a/2
c     Vy(2) = 0.0
c     Vx(3) = -1*a/2
c     Vy(3) = -1*V0
c     Vx(4) = a/2
c     Vy(4) = -1*V0
c     Vx(5) = a/2
c     Vy(5) = 0.0
c     Vx(6) = 20.0
c     Vy(6) = 0.0
c     call pgline (6,Vx,Vy)

      x=0.0
      Vx(1)=-1*V0
      call pgpt(1,x,Vy(1),-1)
      do i = 1, n_steps
        x = i*delta_x
        Vy(1) = -1*V0 + x*x
        call pgpt(1,x,Vy(1),-1)
        x = -1*x
        call pgpt(1,x,Vy(1),-1)
      end do

      i_colour = 2
 
 100  write(*,*) 'Type in parity (1=Odd;0=Even) and trial eigenvalue:'
      read (*,*) odd_even, eigenvalue
c     -38.7315 for example

      if ((odd_even.ne.1).and.(odd_even.ne.0)) goto 100

      if (odd_even .eq. 1) then
        phi_ini = 0.0
        dphi_dx_ini = 1.0
      end if

      if (odd_even .eq. 0) then
        phi_ini = 1.0
        dphi_dx_ini = 0.0
      end if

      call pgsci(i_colour)
      call pgpt(1,0.0,phi_ini + eigenvalue,-1)
  
      n_steps = int(40.0/delta_x)

      phi = phi_ini
      dphi_dx = dphi_dx_ini
      half_a = a/2

      do i = 1, n_steps

        x = i*delta_x

c       if (x.le.half_a) then
c         V = -1*V0
c        else
c         V = 0.0
c       end if

        V = -1*V0 + x*x

        phi_new = phi + dphi_dx*delta_x
        dphi_dx = dphi_dx + (V - eigenvalue)*phi_new*delta_x
        phi = phi + dphi_dx*delta_x

        call pgpt (1,x,phi + eigenvalue,-1)
        if (odd_even.eq.0) then
          call pgpt(1,-1*x,phi + eigenvalue,-1)
         else
          call pgpt(1,-1*x,-1*phi + eigenvalue,-1)
        end if

      end do

 200  write(*,*) 'Try one more time ? (1=Yes;0=No)'
      read (*,*) i_again
      if (i_again.eq.1) then
        i_colour = i_colour + 1
        goto 100
      end if
      if (i_again.ne.0) goto 200

      call pgclos

      end