PROGRAM eigen
      IMPLICIT NONE
      REAL*8 V0, a, dx, xmax
      INTEGER IR_
      CALL GWopen(IR_, 0)
      CALL parameters(V0,a,xmax,dx)
      CALL plot_potential(REAL(V0),REAL(a),REAL(xmax))
      CALL Euler(V0,a,dx,xmax)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE parameters(V0,a,xmax,dx)
      IMPLICIT NONE
      REAL*8 V0, a, xmax, dx
      WRITE(*,'(A,$)') 'magnitude of well depth = '
      READ(*,*) V0
      WRITE(*,'(A,$)') 'half width of well = '
      READ(*,*) a
      WRITE(*,'(A,$)') 'step size = '
      READ(*,*) dx
      WRITE(*,'(A,$)') 'maximum value of x to be plotted = '
      READ(*,*) xmax
      END
C
      SUBROUTINE plot_potential(V0,a,xmax)
*     井戸型ポテンシャルの描画
      IMPLICIT NONE
      REAL V0, a, xmax, GWaspect
      INTEGER IC1_, IR_
      DATA IC1_/1/
      CALL GWvport(IR_, 0.0, 0.0, GWaspect(1), 0.4)
      CALL GWsetpen(IR_, 13, -1, -1, -1)
      CALL GWindow(IR_, -1.05*xmax, -1.1*V0, 1.05*xmax, 1.1*V0)
      CALL GWsavevp(IR_, IC1_)
      CALL GWsetogn(IR_, 0, IC1_)
      CALL GWline(IR_, -xmax, 0.0, -a, 0.0) ! 水平の直線
      CALL GWline(IR_, -a, 0.0, -a, -V0)    ! 垂直の直線
      CALL GWmove2(IR_, -a, -V0)
      CALL GWline2(IR_, a, -V0)
      CALL GWline2(IR_, a, 0.0)
      CALL GWline2(IR_, xmax, 0.0)
      END
C
      SUBROUTINE Euler(V0,a,dx,xmax)
      IMPLICIT NONE
      REAL*8 E, d2phi, dphi, phi, phi_old, x, x_old, V0, a, dx, xmax, V
      REAL GWaspect
      INTEGER IC2_, IR_, parity
      LOGICAL Ltmp1_, Ltmp2_
      DATA IC2_/2/
      CALL GWvport(IR_, 0.0, 0.4, GWaspect(1), 1.0)
      CALL GWindow(IR_, REAL(-1.01*xmax), -4.0, REAL(1.01*xmax), 4.0)
      CALL GWsavevp(IR_, IC2_)
      CALL GWsetogn(IR_, 0, IC2_)
      WRITE(*,'(A,$)') 'even or odd parity (1 or -1) = '
      READ(*,*) parity
      CALL GWerase(IR_, IC2_, -1)
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         WRITE(*,'(A,$)') 'E = '
         READ(*,*) E
         IF(E .EQ. 0) RETURN
         CALL GWsetpen(IR_, 16, -1, -1, -1)
         IF(parity .EQ. -1) THEN
            phi = 0                         ! x = 0 での初期値
            dphi = 1                        ! 一次導関数
         ELSE
            phi = 1
            dphi = 0
         END IF
         x = 0
         Ltmp2_ = .TRUE.
         DO WHILE(Ltmp2_)
            x_old = x
            phi_old = phi
            x = x + dx
            d2phi = 2*(V(x,V0,a) - E)*phi   ! 無次元の単位系
            dphi = dphi + d2phi*dx          ! オイラー-クロマー・アルゴリズム
            phi = phi + dphi*dx
*           波動関数をプロットする
            CALL GWline(IR_, REAL(x_old), REAL(phi_old), REAL(x)
     +      , REAL(phi))
            CALL GWline(IR_, REAL(-x_old), REAL(phi_old*parity)
     +      , REAL(-x), REAL(phi*parity))
            Ltmp2_ = x .LE. xmax
         END DO
      END DO
      END
C
      FUNCTION V(x,V0,a)
      IMPLICIT NONE
      REAL*8 V, x, V0, a
      IF(abs(x) .GT. a) THEN
         V = V0
      ELSE
         V = 0
      END IF
      END