bvp4c show Singular Jacobian Error Message

Why Program Error -

Error using bvp4c (line 251)

Unable to solve the collocation equations -- a

singular Jacobian encountered.

Error in test_si (line 32)

sol1 = bvp4c(@addictive_ode,@addictive_bc,solinit);

close all;clear all;
% global N mu d u1 u2 k1 k2 k3 bet1 bet2 gam1 gam2 epsil1 epsil2 si;
xa=0;
xb=36;
x = linspace(xa,xb);
si = zeros(size(x));
for i = 1:length(x)
if x(i) >= xa & x(i)<=12;
    si(i)=0.1526.*x(i).^2-2.5131.*x(i)+83.543;
elseif x(i)>12 & x(i)<=24;
    si(i)=0.0005*x(i).^6-0.0171*x(i).^5+0.2288*x(i).^4-1.4429*x(i).^3+4.5116*x(i).^2-7.4376*x(i)+79.707;
else x>24;
    si(i)=-0.0054*x(i).^4+0.1722*x(i).^3-1.7406*x(i).^2+6.3251*x(i)+69.433;
end
end
% si = 0.1526*x.^2-2.5131*x+83.543;
% si = 0.0005*x.^6-0.0171*x.^5+0.2288*x.^4-1.4429*x.^3+4.5116*x.^2-7.4376*x+79.707;
mu = 0.000833;
d = 0.000666;
epsil1 = 0.0020;
epsil2 = 0.000634;
bet1 = 0.002453;
bet2 = si*(0.25*0.02);
gam1 = 0.0048;
gam2 = si*(0.25*0.02)+0.00013;
k1 = 1;
k2 = 1.5;
k3 = 0.01;
% xa=0;
% xb=24;
solinit = bvpinit(linspace(xa,xb,100),[50 60 100 10 15 50]);
sol1 = bvp4c(@addictive_ode,@addictive_bc,solinit);
sol2 = bvp4c(@addictive_ode1,@addictive_bc,solinit);
% x = linspace(xa,xb);
y = deval(sol1,x);
y1 = deval(sol2,x);
for i=1:size(sol2.x,2)
    N=sol2.y(1,i)+sol2.y(2,i)+sol2.y(3,i);
    u1(i)=min(max(0,(((sol2.y(4,i)-sol2.y(5,i))*(sol2.y(1,i)*sol2.y(3,i)))/(2*k2*N))),1.0);
    u2(i)=min(max(0,(((sol2.y(5,i)-sol2.y(6,i))*(sol2.y(2,i)*sol2.y(3,i)))/(2*k3*N))),1.0);
end
figure()
plot(x,y(3,:),'LineWidth',2.5) ,hold on
plot(x,y1(3,:),'-.','LineWidth',2.5)
plot(x,y1(2,:),'-.','LineWidth',2.5)
plot(x,y(2,:),'LineWidth',2.5)
xlabel('time','interpreter','latex')
% title('Evolution of $S$ and $A$','interpreter','latex')
legend('$A$ w/o c.','$A$ with c.','$S$ with c.','$S$ w/o c.',...
    'Location','Best','interpreter','latex')
axis([0 xb 0 1])
figure()
plot(x,y(1,:),'LineWidth',2.5) ,hold on
plot(x,y1(1,:),'-.','LineWidth',2.5)
axis([0 xb 0 .5])
xlabel('time','interpreter','latex')
% title('Evolution of $N$','interpreter','latex')
legend('$N$ w/o c.','$N$ with c.','Location','Best',...
    'interpreter','latex')
% axis([0 24 0 1])
figure()
plot(sol2.x,u1,'LineWidth',2.5)
xlabel('time','interpreter','latex')
% title('Variation of the optimal control $u_1$',...
%     'interpreter','latex')
% axis([0 xb 0 6*10^-3])
figure()
plot(sol2.x,u2,'LineWidth',2.5)
xlabel('time','interpreter','latex')
% title('Variation of the optimal control $u_2$',...
%     'interpreter','latex')
axis([0 xb 0 1])
figure()
plot(x,y1(4,:),'LineWidth',2.5),hold on
plot(x,y1(5,:),'LineWidth',2.5)
plot(x,y1(6,:),'LineWidth',2.5)
legend('p_1','p_2','p_3')
% -------------------------------------------------------------
% -------------------------- Function -------------------------
% -------------------------------------------------------------

function dydx = addictive_ode(x,y)
global N mu d u1 u2 k1 k2 k3 bet1 bet2 gam1 gam2 epsil1 epsil2 si;
N = y(1)+y(2)+y(3);
xa=0;
xb=36;
% si = 0.1526*x.^2-2.5131*x+83.543;
% si = 0.0005*x.^6-0.0171*x.^5+0.2288*x.^4-1.4429*x.^3+4.5116*x.^2-7.4376*x+79.707;
si = zeros(size(x));
for i = 1:length(x)
if x(i) >= xa & x(i)<=12;
    si(i)=0.1526.*x(i).^2-2.5131.*x(i)+83.543;
elseif x(i)>12 & x(i)<=24;
    si(i)=0.0005*x(i).^6-0.0171*x(i).^5+0.2288*x(i).^4-1.4429*x(i).^3+4.5116*x(i).^2-7.4376*x(i)+79.707;
else x>24;
    si(i)=-0.0054*x(i).^4+0.1722*x(i).^3-1.7406*x(i).^2+6.3251*x(i)+69.433;
end
end
mu = 0.000833;
d = 0.000666;
epsil1 = 0.0020;
epsil2 = 0.000634;
bet1 = 0.002453;
bet2 = si*(0.25*0.02);
gam1 = 0.0048;
gam2 = si*(0.25*0.02)+0.00013;
k1 = 1;
k2 = 1.5;
k3 = 0.01;
u1=0;
u2=0;
dydx = [(mu*N)-(d*y(1))-((bet1+u1)*(y(1)*y(3))/N)-(bet2*y(1))+(epsil2*y(3))
    ((bet1+u1)*(y(1)*y(3))/N)+(bet2*y(1))-(d*y(2))-((gam1+u2)*(y(2)*y(3))/N)-(gam2*y(2))+(epsil1*y(3))
    ((gam1+u2)*(y(2)*y(3))/N)+(gam2*y(2))-(epsil2*y(3))-(epsil1*y(3))-(d*y(3))
    -y(4)*(mu-d-(bet1+u1)*(y(3)/N)+(bet1+u1)*(y(1)*y(2)/N^2)-bet2)-...
    y(5)*((bet1+u1)*(y(3)/N)-(bet1+u1)*(y(1)*y(2)/N^2)+bet2+(gam1+u2)*(y(2)*y(3))/N^2)+...
    y(6)*((gam1+u2)*(y(2)*y(3))/N^2)
    -k1-y(4)*(mu+(bet1+u1)*(y(1)*y(3)/N^2))-...
    y(5)*(-(bet1+u1)*(y(1)*y(3)/N^2)-d-(gam1+u2)*(y(3)/N)+(gam1+u2)*(y(2)*y(3)/N^2)-gam2)-...
    y(6)*((gam1+u2)*(y(3)/N)-(gam1+u2)*(y(2)*y(3)/N^2)+gam2)
    -y(4)*(mu-(bet1+u1)*(y(1)/N)+(bet1+u1)*(y(1)*y(2)/N^2)+epsil2)-...
    y(5)*((bet1+u1)*(y(1)/N)-(bet1+u1)*(y(1)*y(2)/N^2)-(gam1+u2)*(y(2)/N)+(gam1+u2)*(y(2)*y(3)/N^2)+epsil1)-...
    y(6)*((gam1+u2)*(y(2)/N)-(gam1+u2)*(y(2)*y(3)/N^2)-epsil2-d-epsil1)];
end

% -------------------------------------------------------------
% ---------------------- Function control ---------------------
% -------------------------------------------------------------

function dydx = addictive_ode1(x,y)
global N mu d u1 u2 k1 k2 k3 bet1 bet2 gam1 gam2 epsil1 epsil2 si;
N = y(1)+y(2)+y(3);
xa=0;
xb=36;
% si = 0.1526*x.^2-2.5131*x+83.543;
% si = 0.0005*x.^6-0.0171*x.^5+0.2288*x.^4-1.4429*x.^3+4.5116*x.^2-7.4376*x+79.707;
si = zeros(size(x));
for i = 1:length(x)
if x(i) >= xa & x(i)<=12;
    si(i)=0.1526.*x(i).^2-2.5131.*x(i)+83.543;
elseif x(i)>12 & x(i)<=24;
    si(i)=0.0005*x(i).^6-0.0171*x(i).^5+0.2288*x(i).^4-1.4429*x(i).^3+4.5116*x(i).^2-7.4376*x(i)+79.707;
else x>24;
    si(i)=-0.0054*x(i).^4+0.1722*x(i).^3-1.7406*x(i).^2+6.3251*x(i)+69.433;
end
end
mu = 0.000833;
d = 0.000666;
epsil1 = 0.0020;
epsil2 = 0.000634;
bet1 = 0.002453;
bet2 = si*(0.25*0.02);
gam1 = 0.0048;
gam2 = si*(0.25*0.02)+0.00013;
k1 = 1;
k2 = 1.5;
k3 = 0.01;
u1=min(max(0,((y(4)-y(5))*y(1)*y(3))/(2*k2*N)),1.0);
u2=min(max(0,((y(5)-y(6))*y(2)*y(3))/(2*k3*N)),1.0);
dydx = [(mu*N)-(d*y(1))-((bet1+u1)*(y(1)*y(3))/N)-(bet2*y(1))+(epsil2*y(3))
    ((bet1+u1)*(y(1)*y(3))/N)+(bet2*y(1))-(d*y(2))-((gam1+u2)*(y(2)*y(3))/N)-(gam2*y(2))+(epsil1*y(3))
    ((gam1+u2)*(y(2)*y(3))/N)+(gam2*y(2))-(epsil2*y(3))-(epsil1*y(3))-(d*y(3))
    -y(4)*(mu-d-(bet1+u1)*(y(3)/N)+(bet1+u1)*(y(1)*y(2)/N^2)-bet2)-...
    y(5)*((bet1+u1)*(y(3)/N)-(bet1+u1)*(y(1)*y(2)/N^2)+bet2+(gam1+u2)*(y(2)*y(3))/N^2)+...
    y(6)*((gam1+u2)*(y(2)*y(3))/N^2)
    -k1-y(4)*(mu+(bet1+u1)*(y(1)*y(3)/N^2))-...
    y(5)*(-(bet1+u1)*(y(1)*y(3)/N^2)-d-(gam1+u2)*(y(3)/N)+(gam1+u2)*(y(2)*y(3)/N^2)-gam2)-...
    y(6)*((gam1+u2)*(y(3)/N)-(gam1+u2)*(y(2)*y(3)/N^2)+gam2)
    -y(4)*(mu-(bet1+u1)*(y(1)/N)+(bet1+u1)*(y(1)*y(2)/N^2)+epsil2)-...
    y(5)*((bet1+u1)*(y(1)/N)-(bet1+u1)*(y(1)*y(2)/N^2)-(gam1+u2)*(y(2)/N)+(gam1+u2)*(y(2)*y(3)/N^2)+epsil1)-...
    y(6)*((gam1+u2)*(y(2)/N)-(gam1+u2)*(y(2)*y(3)/N^2)-epsil2-d-epsil1)];
 end

% -------------------------------------------------------------
% -------------------------- Boundary -------------------------
% -------------------------------------------------------------

function res = addictive_bc(ya,yb)
res = [ya(1)-0.4897
    ya(2)-0.4018
    ya(3)-0.1085
    yb(4)-0
    yb(5)-0
    yb(6)-0];
 end