Matlab bvp5c求解含特征值边值问题时特征值不一致的求助

我正在参考一篇论文，用Matlab的bvp5c函数求解含特征值的边值问题。控制方程为二阶常微分方程，其中α、At为常数，Q(σ)是已知的误差函数，ϕ为未知函数，Ψ是待求的特征值；边界条件设置为左侧固定函数值与导数值，右侧需满足渐近衰减条件（有时需调整右侧边界条件）。

多数情况下代码运行结果和论文图1趋势相似，但渐近值和论文结果不符。最关键的问题是：修改特征值初始猜测值Psi1时，求解结果差异极大——比如初始猜1时，渐近值约11；改成11时，渐近值变成13左右；设为0时，渐近值仅8.5。

我试过用更熟悉的Mathematica求解，但它始终返回Psi的线性解，得不到渐近行为。下面是我的Matlab代码，恳请各位提供解决建议或替代方法：

eigenvalues0 = zeros(1, 30);
global alpha;
global At;
options = bvpset('RelTol',1e-5,'AbsTol',1e-5);
At = 0;
for i=1:30    
alpha=1/i;
Psi1 = 1;
solinit = bvpinit(linspace(-3,6,30000),@mat4init,Psi1);
sol = bvp5c(@mat4ode, @mat4bc, solinit,options);
 eigenvalues0(1,i) = sol.parameters;
end
figure(Color="white");
plot(eigenvalues0);
grid on;
legend('0','0.243','0.5','0.843','0.999','FontSize', 14)
xlabel('\alpha^{-1}','FontSize', 14);
ylabel('\Psi','FontSize', 14);
set(gca, 'FontSize', 14);
function dydx = mat4ode(x,y,Psi1) % equation being solved
global alpha;
global At;
dydx = [y(2)
        (-alpha^2*At*2/sqrt(pi)*exp(-x.^2)*y(2)+(1+At*erf(x)-alpha*Psi1*2/sqrt(pi)*exp(-x.^2))*y(1))/(alpha^2*(1+At*erf(x)))];
end
function res = mat4bc(ya,yb,Psi1) % boundary conditions
global alpha;
res = [ya(2)-1/alpha*exp(-3/alpha)
       abs(yb(1))+abs(yb(2))
       %yb(2)
       ya(1)-exp(-3/alpha)];
end
function yinit = mat4init(x) % initial guess function
global alpha;
yinit = [exp(-abs(x)/alpha)
        -sign(x)*1/alpha*exp(-abs(x)/alpha)];
end

内容的提问来源于stack exchange，提问作者Repentanze