エルミート多項式 [FORTRAN]

こんばんは　管理人のさくらえびです。
今日はエルミート多項式に挑戦。
エルミート多項式と言えば、調和振動子関数でもおなじみ(物理分野だとお馴染みのです。)
さっそく、MathJaxでスタイリッシュに書きましょう。
エルミート多項式は、以下の微分方程式を満たす多項式\(H_n(x)\)です。
\[
(\frac{d^2}{dx^2}-2x\frac{d}{dx}+2n)H_n(x) = 0
\]
決まった！鮮やかに決まった！！
ちなみに、TeX表記だといまのはこんな感じです。

(\frac{d^2}{dx^2}-2x\frac{d}{dx}+2n)H_n(x) = 0

さて、この微分方程式を解くと、下記のようなnごとの式表が出ます。
\(H_0(x)= 1\)
\(H_1(x)= 2x\)
\(H_2(x)= 4x^2-2\)
\(H_3(x)= 8x^3-12x\)
漸化式は、
\[
H_{n+1}(x) = 2xH_n(x) -2nH_{n-1}(x)
\]
なので、そろそろソースコード
fortran

function Hnx(n,x)
implicit none
double precision Hnx,x
double precision tmpH,tmp1,tmp2
integer*4 n,i
if (n == 0) then
Hnx = 1.d0
else if(n == 1) then
Hnx = 2.d0*x
else
tmpH = 1.d0
tmp1 = 2.d0*x
do i = 1,n-1,1
tmp2 = tmp1
tmp1 = 2.d0*x*tmp1 - 2.d0*dble(i)*tmpH
tmpH = tmp2
enddo
Hnx = tmp1
endif
end function Hnx

結果は同じなので、最後にノード2までのグラフを載せておきます。
ちなみに、グラフ外に飛んでいくのが好きじゃないので、\(exp(-x^2/2)\)を掛けています。