Warning: Results may be inaccurate because of an ill-conditioned Jacobian whose reciprocal condition number is 2.03573e-39

I’m trying to solve an boundary value problem with eigenvalues, following the article ( https://doi.org/10.1063/5.0155331 )

where and At are certain numbers and Q() is a known function, erf function in this case.
is the unknown function and the eigenvalue.
The boundary condition is , so I specified
For small alphas my code with bvp5c works well.
However for large alphas, my solution deviates from the result in the article, and I’m encountering messages like
Warning: Results may be inaccurate because of an ill-conditioned Jacobian whose reciprocal condition number is 2.03573e-39
I guess the boundary values may be too small and tried to scale my boundary values with C*exp(-3/alpha) but the problem remains.

How can I solve the problem? Fig.1 in https://doi.org/10.1063/5.0155331 can be referred as a correct solution to this problem. My solution looks similar to it but deviates in asymptotic values.
eigenvalues0999 = zeros(1, 30);

global alpha;
global At;
options = bvpset(‘RelTol’,1e-5,’AbsTol’,1e-5);

At = 0.999;
for i=1:30
alpha=1/i;
Psi1 = 1;
solinit = bvpinit(linspace(-3,6,30000),@mat4init,Psi1);
sol = bvp5c(@mat4ode, @mat4bc, solinit,options);
eigenvalues0999(1,i) = sol.parameters;
end

figure(Color="white");
plot(eigenvalues0999);
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)];
endI’m trying to solve an boundary value problem with eigenvalues, following the article ( https://doi.org/10.1063/5.0155331 )

global alpha;
global At;
options = bvpset(‘RelTol’,1e-5,’AbsTol’,1e-5);

At = 0.999;
for i=1:30
alpha=1/i;
Psi1 = 1;
solinit = bvpinit(linspace(-3,6,30000),@mat4init,Psi1);
sol = bvp5c(@mat4ode, @mat4bc, solinit,options);
eigenvalues0999(1,i) = sol.parameters;
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 I’m trying to solve an boundary value problem with eigenvalues, following the article ( https://doi.org/10.1063/5.0155331 )

global alpha;
global At;
options = bvpset(‘RelTol’,1e-5,’AbsTol’,1e-5);

At = 0.999;
for i=1:30
alpha=1/i;
Psi1 = 1;
solinit = bvpinit(linspace(-3,6,30000),@mat4init,Psi1);
sol = bvp5c(@mat4ode, @mat4bc, solinit,options);
eigenvalues0999(1,i) = sol.parameters;
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 bvp5c, solve MATLAB Answers — New Questions

Cart

Cart