why this error? >> Untitledyasir Error using bvp4c (line 251) Unable to solve the collocation equations -- a singular Jacobian encountered. Error in Untitledyasir (line 55) s

function slipflow

format long g

%Define all parameters

% Boundary layer thickness & stepsize

etaMin = 0;

etaMax1 = 15;

etaMax2 = 15; %15, 10

stepsize1 = etaMax1;

stepsize2 = etaMax2;

% Input for the parameters

A=1; %velocity slip

B=0.2; %thermal slip

S=2.4; %suction(2.3,2.4,2.5)

Pr=6.2; %prandtl number

lambda=-1; %stretching shrinking

a=0.01; %phil-1st nanoparticle concentration

b=0.01; %(0.01,0.05)phi2-2nd nanoparticle concentration

c=a+b; %phi-hnf concentration of hybrid nanoparticle

%%%%%%%%%%% 1st nanoparticle properties (Al2O3)%%%%%%%%%%%%

C1=765;

P1=3970;

K1=40;

B1=0.85/((10)^5);

s1=35*(10)^6; %MHD

%%%%%%%%%%% 2nd nanoparticle properties (Cu)%%%%%%%%%%%%

C2=385; %specific heat

P2=8933; %density

K2=400; %thermal conductivity

B2=1.67/((10)^5); %thermal expansion

s2=(59.6)*(10)^6; %MHD

%%%%%%%%%%% Base fluid properties %%%%%%%%%%%%

C3=4179; %specific heat

P3=997.1; %density

K3=0.613; %thermal conductivity

B3=21/((10)^5); %thermal expansion

s3=0.05; %MHD

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%multiplier%%%%%%%%%%%%%%%%%%%

H1=P1*C1; %pho*cp nanoparticle 1

H2=P2*C2; %pho*cp nanoparticle 2

H3=P3*C3; %pho*cp base fluid

H4=a*H1+b*H2+(1-c)*H3; %pho*cp hybrid nanofluid

H5=a*P1+b*P2+(1-c)*P3; %pho hybrid nanofluid

H6=1/((1-c)^2.5); % mu hybrid nanofluid / mu base fluid

H7=a*(P1*B1)+b*(P2*B2)+(1-c)*(P3*B3); % thermal expansion of hybrid nanofluid

%Kn=K3*(K1+2*K3-2*a*(K3-K1))/(K1+2*K3+a*(K3-K1)); %thermal conductivity of nanofluid

Kh=(((a*K1+b*K2)/c)+2*K3+2*(a*K1+b*K2)-2*c*K3)/(((a*K1+b*K2)/c)+2*K3-(a*K1+b*K2)-2*c*K3); %khnf/kf

H8=(((a*s1+b*s2)/c)+2*s3+2*(a*s1+b*s2)-2*c*s3)/(((a*s1+b*s2)/c)+2*s3-(a*s1+b*s2)-2*c*s3); % \sigma hnf/ \sigma f

D1=(H5/P3)/H6;

D3=(H7/(P3*B3))/(H5/P3); % multiplier of boundary parameter

D2= Pr*((H4/H3)/Kh);

D4=H8/(H5/P3); %multiplier MHD

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% First solution %%%%%%%%%%%%%%%%%%%

options = bvpset('stats','off','RelTol',1e-10);

solinit = bvpinit (linspace (etaMin, etaMax1, stepsize1),@(x)OdeInit1(x,A,S,lambda));

sol = bvp4c(@(x,y)OdeBVP(x,y,Pr,D1,Kh,H4,H3,S), @(ya,yb)OdeBC(ya, yb, A, S, B, lambda), solinit, options);

eta = linspace (etaMin, etaMax1, stepsize1);

y= deval (sol,eta);

figure(1) %velocity profile

plot(sol.x,sol.y(2,:),'-')

xlabel('\eta')

ylabel('f`(\eta)')

hold on

figure(2) %temperature profile

plot(sol.x,sol.y(4,:),'-')

xlabel('\eta')

ylabel('\theta(\eta)')

hold on

% saving the out put in text file for first solution

descris =[sol.x; sol.y];

save 'sliphybrid_upper.txt' descris -ascii

% Displaying the output for first solution

fprintf('\n First solution:\n');

fprintf('f"(0)=%7.9f\n',y(3)); % reduced skin friction

fprintf('-theta(0)=%7.9f\n',-y(5)); %reduced local nusselt number

fprintf('Cfx=%7.9f\n',H6*(y(3))); % skin friction

fprintf('Nux=%7.9f\n',-Kh*y(5)); % local nusselt number

fprintf('\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% second solution %%%%%%%%%%%%%%%%%%%

options = bvpset('stats','off','RelTol',1e-10);

solinit = bvpinit (linspace (etaMin, etaMax2, stepsize2),@(x)OdeInit2(x,A,S,lambda));

sol = bvp4c(@(x,y)OdeBVP(x,y,Pr,D1,Kh,H4,H3,S), @(ya,yb)OdeBC(ya, yb, A, S, B, lambda), solinit, options);

eta= linspace (etaMin, etaMax2, stepsize2);

y = deval (sol,eta);

figure(1) %velocity profile

plot(sol.x,sol.y(2,:),'--')

xlabel('\eta')

ylabel('f`(\eta)')

hold on

figure(2) %temperature profile

plot(sol.x,sol.y(4,:),'--')

xlabel('\eta')

ylabel('\theta(\eta)')

hold on

% saving the out put in text file for second solution

descris=[sol.x; sol.y];

save 'sliphybrid_lower.txt'descris -ascii

% Displaying the output for first solution

fprintf('\nSecond solution:\n');

fprintf('f"(0)=%7.9f\n',y(3)); % reduced skin friction

fprintf('-theta(0)=%7.9f\n',-y(5)); %reduced local nusselt number

fprintf('Cfx=%7.9f\n',H6*(y(3))); % skin friction

fprintf('Nux=%7.9f\n',-Kh*y(5)); % local nusselt number

fprintf('\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end

% Define the ODE function

function f = OdeBVP(x,y,Pr,D1,Kh,H4,H3,S)

f =[y(2);y(3);(1/3)*D1*((y(2)^2)-2*(y(1)*y(3))+2*S);y(5);-2/3*(Pr*((H4/H3)/Kh)*y(1)*y(5))];

end

% Define the boundary conditions

function res = OdeBC(ya, yb, A, S, B, lambda)

res= [ya(1)-S;ya(2)-lambda;ya(4)-1;yb(3)-1;yb(4)];

end

% setting the initial guess for first solution

function v = OdeInit1(x,A,S,lambda)

v=[S+0.56;0;0;0;0];

end

% setting the initial guess for second solution

function v1 =OdeInit2(x, A, S,lambda)

v1 = [exp(-x);exp(-x);-exp(-x);-exp(-x);-exp(-x)];

end