How can I solve the error "Unable to solve the collocation equations -- a singular Jacobian encountered" in bvp4c?

I tried to solve the bvp using bvp4c procedure. But I have an error "Unable to solve the collocation equations -- a singular Jacobian encountered." How to resolve this issue?

proj()

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

Error in solution>proj (line 31)
sol= bvp4c(@projfun,@projbc,solinit,options);

function sol= proj



A = 0.5;
    Gr = 0.7;
    Gc = 0.5;
    Kp = 3.0;
    beta = 0.5;
    Pr = 0.3;
    Df = 0.2;
    Sc = 0.1;
    L0 = 0.5;
    Sr = 0.3;
    M = 1;

    Bi1 = 0.5;
    Bi2 = 0.5;
    K0 = 0.3;

myLegend1 = {};myLegend2 = {};
rr = [0.3 0.5 0.8];
for i =1:numel(rr)
    Re = rr(i);

   y0 = [1, 0, 1, 1, 0, 1, 0];
options =bvpset('stats','on','RelTol',1e-4);
    x = linspace(0,10,500);
   solinit = bvpinit(x,y0);
      sol= bvp4c(@projfun,@projbc,solinit,options);
      disp((sol.y(1,20)))
    figure(1)
    plot(sol.x,(sol.y(6,:)))
    grid on,hold on
    myLegend1{i}=['Pr = ',num2str(rr(i))];
xlabel('eta');
    ylabel('(thetas-thetaf)/thetas');
    i=i+1;
end
figure(1)
legend(myLegend1)
hold on

    function dy= projfun(~,y)
        dy= zeros(7,1);
          dy(1) = y(2);                                                       
          dy(2) = y(3);                                                       
          dy(3) = (0.5*A*y(3) - Gr*y(4) - Gc*y(6) + (M + Kp + A)*y(2)) / (1 + (1/beta));
          dy(4) = y(5);
          dy(5) = (Pr*(0.5*A*y(5)) - Pr*Df*Sc*(0.5*A*y(7) + 2*A*y(6) - L0*y(6))) / (1 - (Pr*Df*Sc*Sr));
          dy(6) = y(7);
          dy(7) = (Sc*(0.5*A*y(7) + 2*A*dy(6)) - Sr*Pr*(0.5*A*y(5) + 2*A*y(4))) / (1 - (Sr*Df*Pr));
    end

function res= projbc(ya,yb)

res= [
        ya(2) - (1 + K0*ya(3)); 
        ya(5) + Bi1*(1 - ya(4)); 
        ya(7) + Bi2*(1 - ya(6));  
        yb(2); 
        yb(4); 
        yb(6);
        yb(7)];
end
end

[Torsten]

How to resolve this issue?

It's not a technical issue. Usually, there is something wrong with the equations or the boundary conditions. So you should compare both of them with the mathematical description of the problem.

Answer 1

[Torsten]

Your ODE system is linear - thus you can determine its general solution having 7 free parameters. But incorporating of your boundary conditions leads to a linear system of equations where the coefficient matrix A is rank-deficient (rank 6 instead of rank 7). Consequently, your system is not solvable.

A = 0.5;
Gr = 0.7;
Gc = 0.5;
Kp = 3.0;
beta = 0.5;
Pr = 0.3;
Df = 0.2;
Sc = 0.1;
L0 = 0.5;
Sr = 0.3;
M = 1;
Bi1 = 0.5;
Bi2 = 0.5;
K0 = 0.3;
Mat = zeros(7);
Mat(1,2) = 1;
Mat(2,3) = 1;
Mat(3,2) = (M + Kp + A) / (1 + 1/beta);
Mat(3,3) = 0.5*A/(1 + 1/beta);
Mat(3,4) = -Gr/(1 + 1/beta);
Mat(3,6) = -Gc/(1 + 1/beta);
Mat(4,5) = 1;
Mat(5,5) = Pr*0.5*A/(1 - Pr*Df*Sc*Sr);
Mat(5,6) = - Pr*Df*Sc*(2*A - L0) / (1 - Pr*Df*Sc*Sr);
Mat(5,7) = - Pr*Df*Sc*0.5*A/ (1 - Pr*Df*Sc*Sr);
Mat(6,7) = 1;
Mat(7,4) = - Sr*Pr*2*A / (1 - Sr*Df*Pr);
Mat(7,5) = - Sr*Pr*0.5*A / (1 - Sr*Df*Pr);
Mat(7,7) = Sc*(0.5*A+2*A) / (1 - Sr*Df*Pr);
syms x1(t) x2(t) x3(t) x4(t) x5(t) x6(t) x7(t)
eqns = [diff(x1);diff(x2);diff(x3);diff(x4);diff(x5);diff(x6);diff(x7)]-Mat*[x1;x2;x3;x4;x5;x6;x7]==[0;0;0;0;0;0;0];
sol = dsolve(eqns,'MaxDegree',4)

sol = struct with fields:
    x2: C3*(exp((t*(2976375*((13875*3^(1/2)*342589119285075140699353^(1/2))/1883723499362492536448 - 2321348855894375/1076413428207138592256)^(1/6) - 9802324*((3364850596725*((...
    x1: C1 + C2*(exp((t*(2976375*((13875*3^(1/2)*342589119285075140699353^(1/2))/1883723499362492536448 - 2321348855894375/1076413428207138592256)^(1/6) - 9802324*((33648505967...
    x3: C6*exp(-(t*(865^(1/2) - 1))/24) + C7*exp((t*(865^(1/2) + 1))/24) - C2*(exp((t*(2976375*((13875*3^(1/2)*342589119285075140699353^(1/2))/1883723499362492536448 - 23213488...
    x4: C2*(exp((t*(2976375*((13875*3^(1/2)*342589119285075140699353^(1/2))/1883723499362492536448 - 2321348855894375/1076413428207138592256)^(1/6) - 9802324*((3364850596725*((...
    x5: C4*exp((t*(9802324*((3364850596725*((13875*3^(1/2)*342589119285075140699353^(1/2))/1883723499362492536448 - 2321348855894375/1076413428207138592256)^(1/3))/960855558009...
    x6: C3*(exp((t*(2976375*((13875*3^(1/2)*342589119285075140699353^(1/2))/1883723499362492536448 - 2321348855894375/1076413428207138592256)^(1/6) - 9802324*((3364850596725*((...
    x7: C4*exp((t*(9802324*((3364850596725*((13875*3^(1/2)*342589119285075140699353^(1/2))/1883723499362492536448 - 2321348855894375/1076413428207138592256)^(1/3))/960855558009...

eqn1 = subs(sol.x2,t,0)-(1+K0*subs(sol.x3,t,0))==0;
eqn2 = subs(sol.x5,t,0)+Bi1*(1-subs(sol.x4,t,0))==0;
eqn3 = subs(sol.x7,t,0)+Bi2*(1-subs(sol.x6,t,0))==0;
eqn4 = subs(sol.x2,t,10)==0;
eqn5 = subs(sol.x4,t,10)==0;
eqn6 = subs(sol.x6,t,10)==0;
eqn7 = subs(sol.x7,t,10)==0;
[A,b] = equationsToMatrix([eqn1,eqn2,eqn3,eqn4,eqn5,eqn6,eqn7])

A = 

b = 

rank(A)

ans = 6

rank([A,b])

ans = 7

[Torsten]

I suspect the reason is that you don't have any boundary condition involving y(1).