Error using bvp4c Unable to solve the collocation equations -- a singular Jacobian encountered. Error in MHDI_ori (line 54) sol = bvp4c (@OdeBVP, @OdeBC, solinit, options

MHDI_ori()

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

Error in solution>MHDI_ori (line 49)
sol = bvp4c (@OdeBVP, @OdeBC, solinit, options);

function MHDI_ori
format long g
%Define all parameters
global R M S Pr a C1 P1 K1 C2 K2 P2 C3 P3 K3 C4 K4 P4 H1 H2 H3 H4 Kn 
global D1 D2 B1 B2 B3
%Boundary layer thickness & stepsize
etaMin = 0;
etaMax1 = 15;
etaMax2 = 15; %15, 10
stepsize1 = etaMax1;
stepsize2 = etaMax2;
%Input for the parameters
R= 0.1; %Stretching/Shrinking parameter
M=0.2;        %Magnetic field parameter (0,0.1,0.2)
%S=1;      %Suction parameter
Pr=6.2;     %Prandtl number
a=0.1;      %phi1- 1st nanoparticle concentration
%%%%%%%%%%%%%%%%%%%%%%   1st nanoparticle properties (Cu)   %%%%%%%%%%%%%%%
C1=385;  %C=specific heat
P1=8933; %P=rho=density
K1=400;  %K=thermal conductivity
%%%%%%%%%%%%%%%%%%%%%%   2nd nanoparticle properties (Al2O3)   %%%%%%%%%%%%
C2=765;
P2=3970;
K2=40;
%%%%%%%%%%%%%%%%%%%%%%   3rd nanoparticle properties (TiO2)   %%%%%%%%%%%%%%%
C3=686.2;   %C=specific heat
P3=4250;    %P=rho=density
K3=8.9538;  %K=thermal conductivity
%%%%%%%%%%%%%%%%%%%%%%   base fluid   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C4=4179;
P4=997.1;
K4=0.613;
%%%%%%%%%%%%%%%%%%%%%%%%%% multiplier %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H1=P1*C1;               %pho*cp nanoparticle 1
H2=P2*C2;               %pho*cp nanoparticle 2
H3=P3*C3;               %pho*cp nanoparticle 3
H4=P4*C4;               %pho*cp nanoparticle 4
B1=C1/C4;
B2=C2/C4;
B3=C3/C4;
Kn=((K1+2*K4)-2*a*(K4-K1))/((K1+2*K4)+a*(K4-K1)); %thermal conductivity nano
D1=1-a;
D2=D1^(2.5);
%%%%%%%%%%%%%%%%%%%%%%   first solution   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

options = bvpset('stats','off','RelTol',1e-10);     solinit = bvpinit (linspace (etaMin, etaMax1, stepsize1), @OdeInit1);
sol = bvp4c (@OdeBVP, @OdeBC, 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\prime(\eta)')
hold on

figure(2) %temperature profile
plot(sol.x,sol.y(4,:),'-')
xlabel('\eta')
ylabel('\theta(\eta)')
hold on
%Saving the output in txt file for first solution
descris = [sol.x; sol.y];
save 'ainnaz_upper_ori.txt' descris -ascii

%Displaying the output for first solution 
fprintf('\nFirst 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',(1/D1)*y(3));%skin friction
fprintf('Nux = %7.9f\n',-Kn*y(5) );%local Nusselt number
fprintf('\n');
%%%%%%%%%%%%%%%%%%%%%%   second solution   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

options = bvpset('stats','off','RelTol',1e-10); 
solinit = bvpinit (linspace (etaMin, etaMax2, stepsize2), @OdeInit2);
sol = bvp4c (@OdeBVP, @OdeBC, solinit, options);
eta = linspace (etaMin, etaMax2, stepsize2);
y = deval (sol, eta);

figure(1)
plot(sol.x,sol.y(2,:),'--')
xlabel('\eta')
    ylabel('f`(\eta)')
    hold on
    
    figure(2)
    plot(sol.x,sol.y(4,:),'--')
    xlabel('\eta')
    ylabel('\theta(\eta)')
    hold on
    %Saving the output in txt file for first solution
    descris = [sol.x; sol.y];
    save 'ainnaz_lower_ori.txt' descris -ascii
    
    %Displaying the output for second 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', (1/D1)*(y(3)));  %skin friction
    fprintf('Nux = %7.9f\n', -Kn*y(5));       %local Nusselt number
    fprintf('\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end 
%Define the ODE function
function ff = OdeBVP (x, y);%, M, Pr, Kn, S, R, a, D1, D2, B1, H1, H4)
global M Pr Kn S R a D1 D2 B1 H1 H4
ff = [y(2)
      y(3)
      D2*(D1+a*B1)*(-(y(1)*y(3))+y(2)*y(2)-1-M+M*y(2))
      y(5)
      -y(1)*y(5)*Pr*(D1+a*H1/H4)/(Kn)
] ;
end
%Define the boundary condition
function res = OdeBC (ya, yb);%, S, R)
global S R
res = [ya(1)
       yb(2)-R
       ya(4)-1
       yb(2)-1
       yb(4)
       ] ;  
end

%Setting the initial guess for first solution  
function v = OdeInit1 (x);%, S, R)
global S R

v = [0
    -1.25
    0
    0
    0
    ];  
end

%Setting the initial guess for second solution 
function v1= OdeInit2 (x)
v1 = [exp(-x)
    exp(-x)
    exp(-x)
    exp(-x)
    -exp(-x)];
end

Answer 1

[Torsten]

You say you want to have yb(2) - R = 0 and yb(2) - 1 = 0 at the same time as boundary conditions. This is not possible.

[Nurain]

thank you :))))