I am trying to solve a system of ordinary differential equations using Bvp4c but I get a message error that a singular Jacobian encountered, below the code I am working on

test3()

Error using bvp4c
Unable to solve the collocation equations -- a singular Jacobian encountered.

Error in solution>test3 (line 49)
sol = bvp4c(@fsode,@fsbc,solinit,options);

function test3
clear
clc
close all
format long
%% ------------------------------------------------------------------------
% Parameters
%--------------------------------------------------------------------------
NP=100; %the number of meshpoints
eta_infty =15;
%% ------------------------------------------------------------------------
for i=1:5
    A=[2 3 4 5 6];
    gamma=A(i);
    eps=0.001;   %A=0.1*[2 4 6 8 9];       %Medium porosity
    Beta=0.5;    %A=0.1*[1 2 3 4 5];       %Beta casson parameter
    Omega=7;     %A=[4 5 6 7 8];           %cylinder curvature parameter  
    M=6;         %A=[2 3 4 5 6];           %magnetic field parameter
    P=5;         %A=[3 4 5 6 7];           %porous medium parameter
    Gr=0.2;      %A=[0 2 3 4 5];           %themal Grashof number
    Gm=0.3;                                %mass Grashof number
    Alpha=20;                              %curvature angel
    Pr=7.56;     %for water7.56            %prentdl number
    Rd=0;        %A=0.1*[0 1 3 5 7];       %Radiation parameter
    Ec=0.01;                               %Eckert number
    Nt=0.3;      %A=0.1*[1 2 3 4 5];       %thermophoresis parameter     
    Nb=0.2;      %A=0.1*[2 3 4 5 6];       %Brownian motion parameter   
    Sc=5;        %A=[5 7 9 11 13];         %Schmidt number
    %gamma=3;    %A=[2 3 4 5 6];
    xi_1=0.001;                       
    xi_2=0.001;
    QM=0.01;                               %mass flux 
    QH=0.01;                               %heat flux 
    %--------------------------------------------------------------------------
    %lamda=0.001;                          %A=[0.001 0.002 0.003 0.004 0.005];lamda parameter
    %Pe=Re*Pr;                             % peclet number
    %Re=0.001;
    %A1=1/Pr;
    %A2=1/Sc;
    %--------------------------------------------------------------------------
    %  Solve the problem
    %--------------------------------------------------------------------------
    int=ones(1,7); %We have 7 variables
    int=0.*int;
    % int=[0 0 -0.003 -0.35 0 0 -0.05];
    options = [];
    solinit = bvpinit(linspace(0,eta_infty,NP),int);
    sol = bvp4c(@fsode,@fsbc,solinit,options);
    eta = sol.x; %this is eta
    f = sol.y;   %These are the dependenat variables
    FF=f(1,:);
    DF=f(2,:);
    %DDF(i)=-f(3,1)
    TH=f(4,:);
    %DTH(i)=f(5,1);
    Fi =f(6,:);
    %QH=(1+xi_2)*gamma+(1/Pr)*(1+((4*Rd)/3))*DTH;
    %% ------------------------------------------------------------------------
    % Plot the results
    % The model equations
    %--------------------------------------------------------------------------
end
    function dfdeta = fsode(eta,Y)
        F(1) = Y(2);
        F(2) = Y(3);
        F(3) =(-eps*((1+(1/Beta))*(1+(2*Omega*eta))))*(((2*Omega*(1+(1/Beta))*Y(3))/eps)+(1/eps^2)*(Y(1)*Y(3)-(Y(2))^2)-(M+P)*Y(2)+(Gr*Y(4)-Gm*Y(6))*cosd(Alpha));
        F(4)=Y(5);
        F(5)=(-Pr/((1+(2*Omega*eta))*(1+((4/3)*Rd))))*(((2*Omega*(1+((4/3)*Rd))*Y(5))/Pr)+((1+(1/Beta))*(1+(2*Omega*eta))*Ec*(Y(3))^2)+((1+(2*Omega*eta))*((Nt*(Y(5))^2)-(Nb*Y(5)*Y(7))))+(M*Ec*(Y(2))^2)+(Y(1)*Y(5)));
        F(6)=Y(7);
        F(7)=(-1/(1+(2*Omega*eta)))*((2*Omega*Y(7))-((Nt/Nb)*(((1+(2*Omega*eta))*F(5))+(2*Omega*Y(5))))+(Sc*Y(1)*Y(7)));
        dfdeta = [ F(1)
            F(2)
            F(3)
            F(4)
            F(5)
            F(6)
            F(7)
            ];
    end
    % --------------------------------------------------------------------------
    function res = fsbc(ya,yb)
        res = [ya(1)-gamma        % f-gamma=0
            ya(2)                      % f'=0
            ya(4)-1                    %TH=1
            ((Sc/eps)*(ya(6)-xi_1)*ya(1))+ya(7)+((Nt/Nb)*ya(5))-QM       %Mass flux
            yb(2)                      %f'=0
            yb(7)-1                    %phi=1
            (((1+(2*Omega*eta_infty))^(-1/2))*(yb(4)+xi_2)*yb(1))+((1/Pr)*((1+(2*Omega*eta_infty))^(-1/2))*(1+((4/3)*Rd))*yb(5))-QH     %Heat flux
            ];
    end
end

[Torsten]

This error message usually appears if there is something wrong in the setup of your model equations and/or boundary conditions. So we cannot help in this case apart from saying: Recheck your model equations and boundary conditions and their implementation in the code.

[Nora Fattah Ahmed]

Okay.Thank you.