Se for utilizar o material, favor citar como:

1)J.P. Braga e F.S. Carvalho - Métodos Numéricos em Química Quântica, 2021.
2)J.P. Braga, Fundamentos de Química Quântica, Editora UFV, 2007.

clc,clear all,close all
%% colisao H---H
mu=0.503985;
%% Energia
E=0.103675554019775; %%eV
E=E/27.2107; %%au
K2=2*mu*1822.88*E;
%K2=25.0000; %% Angs-2
%% cond inicial e parametros
X(1)=1e-15;h=0.001;
Xmax=15;l=0;
%% Chamando Numerov
[X,R,T,U] = RNM(K2,l,X,h,Xmax,mu);
%% calculando matriz K
[K,Z]=MK(K2,l,X,h,Xmax,T,R,mu);
%% calculando matriz S
[RS,IS]=MS(K);
[RS(length(RS)) IS(length(IS))]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Funções

function [X,R,T,U] = RNM(K2,l,X,h,Xmax,mu)
%% Valores iniciais
i=1;
R(i)=10^35;
T(i)=-((h^2)/12)*Q(K2,l,X(i),mu);
U(i)=(2+10*T(i))/(1-T(i));
%% propagando a solucao de NUMEROV
while X < Xmax
    i=i+1;
    X(i)=X(i-1)+h;
    T(i)=-((h^2)/12)*Q(K2,l,X(i),mu);
    U(i)=(2+10*T(i))/(1-T(i));
    R(i)=U(i) - 1/R(i-1); 
end
end


function Q = Q(k2,l,x,mu)
%% au-2
Q = k2 - l*(l+1)./x.^2 -2*mu*1822.88*pot(x);
%% Angs-2
%Q = k2 - l*(l+1)./(x.^2) - pot(x);
end


function V=pot(x)
%% Calculando pot;
De= 0.173131054131594;
Re=0.74*1.8897259886;
alfa=2.4/1.8897259886;
%% au-2
V=De*((1-exp(-alfa*(x-Re))).^2)-De;
%% Angs-2
%V=1136.0*exp(-2.4*(x-0.74)).*(exp(-2.4*(x-0.74))-2.0);
end


function nl=Bessel1(l,x)
if l==0
    nl=-cos(x);
elseif l==1
    nl=-sin(x) - cos(x)./(x);
else
    nl2(1)=-cos(x);
    nl2(2)=-sin(x) - cos(x)./(x);
    for i=3:l+1
        nl2(i) = (2*(i-2)+1)*nl2(i-1)./(x) - nl2(i-2);
    end
    nl=nl2(i);
end
end


function jl=Bessel2(l,x)
N=100;
J(N+1)=1;
J(N+2)=0;
while N >= 1 
    J(N) = ((2*N+1)/x)*J(N+1) - J(N+2);
    N=N-1;
    
end
p=sin(x)/J(1);
jl=p*J(l+1);
end


function [K,Z]=MK(K2,l,X,h,Xmax,T,R,mu)
Z=sqrt(K2)*X;
% definindo coordenada extra
Zp=sqrt(K2)*(X(length(X))+h);
% acrescentando coordenada extra no vetor de coordenadas
Z=[Z Zp];
% calculando valor de T para coordenada extra, Tp
Tp=-((h^2)/12)*Q(K2,l,X(length(X))+h,mu);
% acrescentando valor extra no vetor T
T=[T Tp];
%% Calculando nl
%B1=Bessel1(l,Z); %% B1=G1
for i=1:size(Z,2)
    B1(i)=Bessel1(l,Z(i)); %% B1=G1
end
%% Calculando jl
for i=1:size(Z,2)
    B2(i)=Bessel2(l,Z(i)); %% B2=F1
end
%% calculando K
I=ones(1,size(Z,2));
JN=(I-T).*B2; %% JN=F1
NN=(I-T).*B1; %% NN=G1
for i=1:size(Z,2)-1
    K(i)=-(R(i)*JN(i)-JN(i+1))./(R(i)*NN(i) - NN(i+1));
end
end


function [RS,IS]=MS(K)
I=ones(1,size(K,2));
RS=(I-K.^2)./(I+K.^2);
IS=-2*K./(I+K.^2); % ou IS=-2K./(I+K.^2);
end