function [DE,NumNode,Psin] = RNNumerov(Qn,dx)

%%
%繰りこみヌメロフ法による1次元時間に依存しないシュレーディンガー方程式の数値解放
%参考：Johnson, Bernard R. "New numerical methods applied to solving the one‐dimensional eigenvalue problem."
%      The Journal of Chemical Physics 67.9 (1977): 4086-4093.
%スクリプト内コメントの式番号はこの論文の式番号
%
%製作者：堀越宗一
%製作日：2019/08/09
%最終更新日：2019/08/13

%%
%Qnのみで与えられる要素を計算

Tn=-(dx^2/12)*Qn;%(32)式
Un=(2+10*Tn)./(1-Tn);%(35)式
An=(0.5-Tn)./(1-Tn);%(44)式

%%
%R,Lの計算と接合点xmの決定

Numx=length(Qn);

Ln=zeros(1,Numx);
Ln(1)=Inf;%x=xminでψ(xmin)=ψ'(xmin)=0;

Rn=zeros(1,Numx);
Rn(Numx)=Inf;%x=xmaxでψ(xmax)=ψ'(xmax)=0;

%x=xmax-1からxMまで内側に積分
for k=Numx-1:-1:1
    Rn(k)=Un(k)-1/Rn(k+1);%(46)式
    if Rn(k)<=1
        m=k;%接合点の決定：%(47)式の3行下に書いてある
        break
    end
end

%左端まで到着した場合適切なxMが見つからなかったのでここで終了
if k==1
    DE=NaN;
    NumNode=NaN;
    Psin=NaN;
else
    
    %適切なxMが見つかった場合、x=xmin+1からxMまで外側に積分
    for k=2:1:m
        Ln(k)=Un(k)-1/Ln(k-1);%(37)式
    end
    
    %%
    %Fnの計算
    
    Fn=zeros(1,Numx);
    Fn(m)=1;%x=xMでFn=1とする（任意の値でよい）：4090ページⅢ章2段落8行目に書いてある
    
    %x=xM-1からx=xminまでFnを計算
    for k=m-1:-1:1
        Fn(k)=Fn(k+1)/Ln(k);%(36)式
    end
    
    %x=xM+1からx=xmaxまでFnを計算
    for k=m+1:1:Numx
        Fn(k)=Fn(k-1)/Rn(k);%(45)式
    end
    
    %%
    %波動関数ψを計算
    
    Psin=Fn./(1-Tn);%(33)式
    Nc=trapz(Psin.^2);%規格化定数
    Psin=Psin/sqrt(Nc);%規格化された波動関数
    
    %%
    %節の数をカウント
    
    Node=diff(sign(Psin));
    NumNode=length(nonzeros(Node(2:m-1)));
    
    %%
    %接合点での整合性のエラー値
    
    DE=(An(m+1)*(1/Rn(m+1)-Ln(m))-An(m-1)*(Rn(m)-1/Ln(m-1)))*(1-Tn(m));%(48)式
    
end

end