MATLAB: Singular Jacobian with BVP4c used to solve eigenvalue problem

Question

I have been trying to solve an eigenvalue problem of the form:

y"+y/x+(Lambda^2)*(1-x^2)y=0

which is a Sturm-Liouville problem that occurs e.g. in the case of heat convection in tubes and channels. The task here is to find the eigenvalues Lambda in order to describe the temperature in the fluid.

As I understand one way to solve this problem is by intorudcing the parameter Lambda into the definition of the differential equation:

function dydx=evp1(x,y,lambda); 
    dydx=[y(2)
         -y(2)/x-lambda*(1-x^2)*y(1)];
end
function res=evp1bc1(ya,yb,lambda); 
    res=[ya(1)-1;yb(1);ya(2)];     
end
function v=evp1Guess1(x); 
    v=[cos(x);sin(x)]; 
end
clear all; close all; echo on; format long;
a=0; b=1; N=10; lambda=0.5;
solInit1=bvpinit(linspace(a,b,N),@evp1Guess1,lambda);
solN1=bvp4c(@evp1,@evp1bc1,solInit1);

This short and elegant code by Zaitsev (in the Book: Handbook of Ordinary DIfferential Equations) finds the eigenvalues for simpler eigenvalue problems (e.g. y"+Lambda.y=0). However, for my problem it creates an error message referring to a singular Jacobian which I assume is due to the presence of the y(1)/x term in the equation at x=0. Is there a way to get around this problem?

Answer 1

Use a=(a small value) instead of a=0.

Best wishes

Torsten.