結果

入力例　1s,2s,2p,3p軌道に固有エネルギーの理論値を代入しています。

縦軸も横軸も無い駄目駄目グラフですが、水素原子の1s,2s,2p,3p軌道の動径波動関数を書いています。
横軸の一目盛りは１ボーア半径、縦軸の単位は任意です。エネルギーの単位はハートリー（２リードベルグ）なので、主量子数nに対してとなります。

s軌道(l=0)は原点で有限の値を持ちますが、p軌道(l=1)は原点で０になります。本来、境界条件から波動関数は遠方では０に収束するはずですが、微分方程式の解法が素朴なオイラー法なので、1s軌道の計算は誤差が溜まって距離が大きいところで発散しています。

MODULE m_plot
USE uhoplot
IMPLICIT NONE
INTEGER, PARAMETER :: nx = 800, ny = 600, noff = 50
CONTAINS
!------------------------------------------------------
SUBROUTINE axis(imax)
IMPLICIT NONE
INTEGER, INTENT(IN) :: imax
INTEGER :: i, ix, iy
CALL gr_move(      noff, ny / 2   ) ! x-axis
CALL gr_line( nx - noff, ny / 2   )
CALL gr_move(      noff,      noff) ! y-axis
CALL gr_line(      noff, ny - noff)
CALL gr_show()
DO i = 1, imax
 ix = (nx - 2 * noff) * i / imax + noff
 iy = ny / 2
 CALL gr_move(ix, iy)
 CALL gr_line(ix, iy + 10)
END DO
RETURN
END SUBROUTINE axis
!------------------------------------------------------
SUBROUTINE line(p1, p2, q1, q2, xmin, xmax, ymin, ymax) 
REAL, INTENT(IN) :: p1, p2, q1, q2, xmin, xmax, ymin, ymax
INTEGER :: ip1, ip2, iq1, iq2
ip1 = (nx - 2 * noff) * p1 / (xmax - xmin) + noff
ip2 = (ny - 2 * noff) * p2 / (ymax - ymin) + ny / 2
iq1 = (nx - 2 * noff) * q1 / (xmax - xmin) + noff
iq2 = (ny - 2 * noff) * q2 / (ymax - ymin) + ny / 2
CALL gr_move(ip1, ip2)
CALL gr_line(iq1, iq2)
RETURN
END SUBROUTINE line
!------------------------------------------------------
END MODULE m_plot
!======================================================
PROGRAM hatom
USE m_plot
IMPLICIT NONE
REAL, PARAMETER :: dr = 0.01, rmin = 0.01, rmax = 20.0
INTEGER :: l 
REAL :: e, r, x, x1, x2
REAL :: a1, b1, a2, b2
l = 0 ! angular momentum
CALL gr_on('Wavefunction of Hydrogen', nx, ny)
CALL axis(INT(rmax))
DO
 WRITE(*, *) ' Input Angular Momentum L & Energy e (e > 100 to exit) '
 READ(*, *) l, e
 IF (e > 100.0) EXIT
! Initial wave function around the origin
 r  = rmin
 x  = r**(l + 1)    ! x = P(r) = rR(r)
 x1 = (l+1) * r**l  ! dx/dr 
 a1 = r
 b1 = x / r
 DO
  CALL schroed( x2, x, e, r, l )
  r  = r + dr
  x1 = x1 + dr * x2
  x  = x  + dr * x1
  a2 = r
  b2 = x / r
  CALL line(a1, b1, a2, b2, 0.0, rmax, 2.0, -2.0)
  a1 = a2
  b1 = b2
  IF (r > rmax) EXIT
 END DO
 CALL gr_show()
END DO
CALL gr_off(0)
STOP
CONTAINS
!------------------------------------------------------
 SUBROUTINE schroed(x2, x, e, r, l)
 REAL   , INTENT(IN ):: x, e, r
 REAL   , INTENT(OUT):: x2
 INTEGER, INTENT(IN) :: l
 REAL :: v, a, b
 v = - 1.0 / r           ! Coulomb potential in Atomic Unit
 a = 2.0 * (v - e)       ! d^2X(r) / dr^2 = ( 2(V(r) - E ) + L(L + 1)/r**2 ) * X(r) 
 b = l * (l + 1) / r**2
 x2 = (a + b) * x 
 RETURN
 END SUBROUTINE schroed
!------------------------------------------------------
END PROGRAM hatom