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
E=70.9; %%eV
Ecol=E/27.2107; %%au
l=0;X=0;
h=0.0001;Xmax=16;
[X,RN,TN,UN] = RNV(Ecol,l,X,h,Xmax);
K=KM(Ecol,l,X,h,TN,RN);
[SR,SI]=SM(K);
S12=SR(1,2)+i*SI(1,2)
S11=SR(1,1)+i*SI(1,1);
S21=SR(2,1)+i*SI(2,1);
S22=SR(2,2)+i*SI(2,2);

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

function [X,RN,TN,UN] = RNV(Ecol,l,X,h,Xmax)
%% R inicial
RN=eye(2)*10^150;
R0=RN^(-1);
%% propagando a solucao de NUMEROV
while X < Xmax
    X=X+h;
    TN=-((h^2)/12)*Q(Ecol,l,X);
    UN=(2*eye(2)+10*TN)*((eye(2)-TN)^(-1));
    RN=UN - R0;
    R0=RN^(-1);    
end


function Q = Q(Ecol,l,x)
%% mu/hbar^2 = 6089.0
U=pot(x);
%% au-2
for i=1:2
    for j=1:2
        Q(i,j) = (2.0*6089.0*Ecol - l*(l+1)./x.^2)*eq(i,j) - U(i,j);
    end
end


function U=pot(x)
U(1,1)=2.0*6089.0*(21.1/x)*exp(-x/0.678);
U(2,2)=2.0*6089.0*(((21.1/x)-12.1)*exp(-x/0.678)+0.617403);
U(1,2)=2.0*6089.0*0.170*exp(-x/0.667);
U(2,1)=U(1,2);


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


function jl=Bessel2(l,x)
Nmax=l;FPA=1;
teste=100;
while teste>1e-6    
N=Nmax;    
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
Nmax=Nmax+20;
p=sin(x)/J(1);
jl=p*J(l+1);
teste=abs(1-abs(FPA)/abs(jl));
FPA=jl;
end


function K=KM(Ecol,l,X,h,TN,RN)
%% determinando valores necessarios no ponto final e um ponto seguinte
E(1,1)=0;E(2,2)=0.617403;
for i=1:2
    Zmax(i,i)=sqrt(2.0*6089.0*(Ecol-E(i,i)))*X;
    Zp(i,i)=sqrt(2.0*6089.0*(Ecol-E(i,i)))*(X+h);
end
Tp=-((h^2)/12)*Q(Ecol,l,X+h);
%% Calculando elementos nl
for i=1:2
    B1(i,i)=Bessel1(l,Zmax(i,i));
    B1p(i,i)=Bessel1(l,Zp(i,i));
end
%% Calculando elementos jl
for i=1:2
    B2(i,i)=Bessel2(l,Zmax(i,i));
    B2p(i,i)=Bessel2(l,Zp(i,i));
end
%% calculando matrizes Nl e Jl
F1=(eye(2)-TN)*B2*(sqrt(Zmax/X))^(-1); %% JN=F1
G1=(eye(2)-TN)*B1*(sqrt(Zmax/X))^(-1); %% NN=G1
F1p=(eye(2)-Tp)*B2p*(sqrt(Zp/(X+h)))^(-1); %% JN=F1
G1p=(eye(2)-Tp)*B1p*(sqrt(Zp/(X+h)))^(-1); %% NN=G1
%% calculando K
K=-((RN*G1 - G1p)^(-1))*(RN*F1-F1p);


function [SR,SI]=SM(K)
SR=((eye(2)+K*K)^(-1))*(eye(2)-K*K);
SI=-2*((eye(2)+K*K)^(-1))*K;