Error solving Boundary value Problem using bv4pc function

Hi, I am getting the following error while solving this BVP.

Error using bvp4c (line 251)
Unable to solve the collocation equations -- a singular Jacobian encountered.

Error in MaxwellYangHWDiff (line 33)
sol = bvp4c(@(r,y)odefun(r,y,phi,chi),@(ycent,ysurf)bcfun(ycent,ysurf,c1,c2,c3,c4),init);

My code is as follows:

clc
clear all
 
R = 0.04; % m Radius of the cylinder
f = 2800 * 10^6; % Hertz Frequency of microwave
c = 2.998 * 10^8; % m/s speed of light
mu_o = 1.257 * 10^-6; % Henry/m free space permeability
eps_o = 8.85 * 10^-12; % Frard/m free space permittivity
k1 = 42.6; % dielectric constant
k2 = 13.1; % dielectric loss 
P_o = 250; % Watts
 
a_o = 2*pi*f/c; %free space wavenumber
E_o = sqrt(P_o*pi*a_o*R/(c*eps_o)); % V/m intensity of incident electric field
phi = R^2*(2*pi*f)^2*mu_o*eps_o*k1; 
chi = R^2*(2*pi*f)^2*mu_o*eps_o*k2;
 
nt = 10; %number of nodes
dR = R/(nt-1); %increment
rstar = 0:(dR/R):(R/R);
 
c1 = R*a_o*((besselj(1,a_o*R)*besselj(0,a_o*R)+bessely(1,a_o*R)*bessely(0,a_o*R))/((besselj(0,a_o*R))^2+(bessely(0,a_o*R))^2));
c2 = 2/(pi*((besselj(0,a_o*R))^2+(bessely(0,a_o*R))^2));
c3 = -(4/pi)*bessely(0,a_o*R)/((besselj(0,a_o*R))^2+(bessely(0,a_o*R))^2);
c4 = -(4/pi)*besselj(0,a_o*R)/((besselj(0,a_o*R))^2+(bessely(0,a_o*R))^2);
 
 
init = bvpinit(linspace(0,0.04,10),[0 0 0 0]);
sol = bvp4c(@(r,y)odefun(r,y,phi,chi),@(ycent,ysurf)bcfun(ycent,ysurf,c1,c2,c3,c4),init);
x = linspace(0,0.04,100); BS=deval(sol,x);
plot(x,BS(1,:))

function dydr = odefun(r,y,phi,chi)
dydr=zeros(4,1);
dydr(2)= -((1/r)*y(2)+phi*y(1)-chi*y(3));
dydr(4)= -((1/r)*y(4)+phi*y(3)-chi*y(1));
end

function res = bcfun(ycent,ysurf,c1,c2,c3,c4)
res =[ycent(2); ycent(4); ysurf(2)+c1*ysurf(1)+c2*ysurf(3)-c3; ysurf(4)+c1*ysurf(3)-c2*ysurf(1)-c4];
end

Answer 1

[darova]

Try r = 1e-3

I think you forgot about dydr(1) and dydr(3)

function dydr = odefun(r,y,phi,chi)
dydr=zeros(4,1);
% dydr(1) = y(2);
dydr(2)= -((1/r)*y(2)+phi*y(1)-chi*y(3));
% dydr(3) = y(4);
dydr(4)= -((1/r)*y(4)+phi*y(3)-chi*y(1));
end

[CM007]

hi, thank you for your answer. I made the changes you suggested. The error I mentioned earlier no longer exists. However, I am getting the following error:

Attempted to access x(1); index out of bounds because numel(x)=0.

Error in bvpinit (line 71)
if x(1) == x(N)

Error in MaxwellYangHWDiff (line 28)
init = bvpinit(linspace(0,0.04,1e-3),[0 0 0 0]);

I tried using different intial points, but it still gives me this error.

Also, there is a correction.

function dydr = odefun(r,y,phi,chi)
dydr=zeros(4,1);
dydr(1) = y(2);
dydr(2)= -((1/r)*y(2)+phi*y(1)-chi*y(3));
dydr(3) = y(4);
dydr(4)= -((1/r)*y(4)+phi*y(3)+chi*y(1));
end

[darova]

I tried this

init = bvpinit(linspace(eps,0.4,100),[0 0 0 0]);
sol = bvp4c(@odefun, @bcfun, init);
plot(sol.x,sol.y)

And get this

See attached script

[CM007]

Thank you! I really appreciate your help!