Anyone who can help on the bvp4c code?

Question

0 votes

The models to be solved are for the Maxwell fluid:

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

Error in solution>Honorsboundary01 (line 46)
sol=bvp4c(@odefun,@BCfun,solinit,options);
function Honorsboundary01
%Special variables
M=1;   %Magnetic parameter
gamma=1;  % Porosity parameter
Rd=1;  %Radiation parameter
lambda=1; %Heat Souce parameter
Pr=1;     %Prandtl number 
Ec=2; %Eckert number 
De=2; %Deboarh number has a range of 0<De<0.1 else the solution may explode
GrT=2;%Thermal Grashof number holds for range 0<GrT<1
GrC=2;%Concentration Grashof number
Le=2;% Lewis number
%System of first order equations
    function dydx=odefun(~,y)
     dydx=zeros(7,1);
     dydx(1)=y(2);
     dydx(2)=y(3);
     dydx(3)=((M+gamma)*y(2)+(y(2)^2)-(M.*De+1)*y(1)*y(3)-GrT*y(4)-GrC*y(6)-2*De*y(1)*y(2)*y(3))/(1-De*(y(1)^2));
     dydx(4)=y(5);
     dydx(5)=(y(4)*y(2)-y(1)*y(5)-Ec*(y(3)^2)-M.*Ec*(y(2)^2)-lambda*y(4))/((1/Pr)*(1+(4/3)*Rd));
     dydx(6)=y(6);
     dydx(7)=Le*y(1)*y(7);
    end
%Boundary conditions of the study
    function Res=BCfun(ya,yb)
      Res=zeros(7,1);
      Res(1)=ya(1);
      Res(2)=ya(2)-1;
      Res(3)=ya(4)-1;
      Res(4)=ya(6)-1;
      Res(5)=yb(2);
      Res(6)=yb(4);
      Res(7)=yb(6);
   
    end
%Initial guess for the problem
solinit=bvpinit(linspace(0,20,10),[0 0 0 0 0 0 0]);
% Vector used for studying the parameter variations
%Implementation of the bvp5c
    options=bvpset('RelTol',1e-10,'AbsTol',1e-10);
sol=bvp4c(@odefun,@BCfun,solinit,options);
eta=linspace(0,10,400);
y=deval(sol,eta);
% 
% Cfx(j)=y(3,1);
% Nux(j)=-(1+(4/3)*Rd(i))*y(5,1);
% Shx(j)=-y(7,1);
    
%  Plots
figure(1);
plot(eta,y(1,:),Linestyle='--',LineWidth=2);
xlabel("\eta");
ylabel("f(\eta)");
hold on
figure(2);
plot(eta,y(2,:),Linestyle='-',LineWidth=2);
xlabel("\eta");
ylabel("f'(\eta)");
hold on
figure(3);
plot(eta,y(3,:),Linestyle='-',LineWidth=2);
xlabel("\eta");
ylabel("f''(\eta)");
hold on
figure(4);
plot(eta,y(4,:),Linestyle='-',LineWidth=2);
xlabel("\eta");
ylabel("\theta(\eta)");
hold on
 figure(5);
 plot(eta,y(5,:),Linestyle='-',LineWidth=2);
 xlabel("\eta");
 ylabel("\theta'(\eta)");
 hold on
%  figure(6)
%  plot(M,Cfx,Cl(i),LineWidth=2);
%  xlabel('M');
%  ylabel("$(1/2)Re_{x}^{(1/2)}$",Interpreter="latex");
%  hold on
%  figure(7)
%  plot(M,Nux,Cl(i),LineWidth=2);
%   xlabel('M');
%  ylabel("$Re_{x}^{-(1/2)}Nu_{x}$",Interpreter="latex");
%  hold on
%  figure(8)
%  plot(M,Shx,Cl(i),LineWidth=2);
%  xlabel('M');
%  ylabel("$Re_{x}^{-(1/2)}Sh_{x}$",Interpreter="latex");
%  hold on
    
figure(1);
legend([repmat('M=',length(M),1),num2str(M')]);
figure(2);
legend([repmat('M=',length(M),1),num2str(M')]);
figure(3);
legend([repmat('M=',length(M),1),num2str(M')]);
figure(4);
legend([repmat('M=',length(M),1),num2str(M')]);
figure(5);
legend([repmat('M=',length(M),1),num2str(M')]);
%fprintf('%.10f\n ',y(5,1));
end

2 Comments
Show None Hide None

Torsten on 26 Apr 2026

Open in MATLAB Online

I didn't test if this solves the problem, but

dydx(6)=y(6);

should be

dydx(6)=y(7);

Walter Roberson on 26 Apr 2026

Open in MATLAB Online

Testing what @Torsten found.

If the NMax parameter of bvpset is configured to more than roughly 2e4, then insetad of mesh tolerance, you get a singular jocobian.

Honorsboundary01()

Warning: Unable to meet the tolerance without using more than 1428 mesh points.
The last mesh of 730 points and the solution are available in the output argument.
The maximum residual is 111.107, while requested accuracy is 1e-10.

function Honorsboundary01

%Special variables

M=1; %Magnetic parameter

gamma=1; % Porosity parameter

Rd=1; %Radiation parameter

lambda=1; %Heat Souce parameter

Pr=1; %Prandtl number

Ec=2; %Eckert number

De=2; %Deboarh number has a range of 0<De<0.1 else the solution may explode

GrT=2;%Thermal Grashof number holds for range 0<GrT<1

GrC=2;%Concentration Grashof number

Le=2;% Lewis number

%System of first order equations

function dydx=odefun(~,y)

dydx=zeros(7,1);

dydx(1)=y(2);

dydx(2)=y(3);

dydx(3)=((M+gamma)*y(2)+(y(2)^2)-(M.*De+1)*y(1)*y(3)-GrT*y(4)-GrC*y(6)-2*De*y(1)*y(2)*y(3))/(1-De*(y(1)^2));

dydx(4)=y(5);

dydx(5)=(y(4)*y(2)-y(1)*y(5)-Ec*(y(3)^2)-M.*Ec*(y(2)^2)-lambda*y(4))/((1/Pr)*(1+(4/3)*Rd));

dydx(6)=y(7);

dydx(7)=Le*y(1)*y(7);

end

%Boundary conditions of the study

function Res=BCfun(ya,yb)

Res=zeros(7,1);

Res(1)=ya(1);

Res(2)=ya(2)-1;

Res(3)=ya(4)-1;

Res(4)=ya(6)-1;

Res(5)=yb(2);

Res(6)=yb(4);

Res(7)=yb(6);

end

%Initial guess for the problem

solinit=bvpinit(linspace(0,20,10),[0 0 0 0 0 0 0]);

% Vector used for studying the parameter variations

%Implementation of the bvp5c

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

sol=bvp4c(@odefun,@BCfun,solinit,options);

eta=linspace(0,10,400);

y=deval(sol,eta);

%

% Cfx(j)=y(3,1);

% Nux(j)=-(1+(4/3)*Rd(i))*y(5,1);

% Shx(j)=-y(7,1);

% Plots

figure(1);

plot(eta,y(1,:),Linestyle='--',LineWidth=2);

xlabel("\eta");

ylabel("f(\eta)");

hold on

figure(2);

plot(eta,y(2,:),Linestyle='-',LineWidth=2);

xlabel("\eta");

ylabel("f'(\eta)");

hold on

figure(3);

plot(eta,y(3,:),Linestyle='-',LineWidth=2);

xlabel("\eta");

ylabel("f''(\eta)");

hold on

figure(4);

plot(eta,y(4,:),Linestyle='-',LineWidth=2);

xlabel("\eta");

ylabel("\theta(\eta)");

hold on

figure(5);

plot(eta,y(5,:),Linestyle='-',LineWidth=2);

xlabel("\eta");

ylabel("\theta'(\eta)");

hold on

% figure(6)

% plot(M,Cfx,Cl(i),LineWidth=2);

% xlabel('M');

% ylabel("$(1/2)Re_{x}^{(1/2)}$",Interpreter="latex");

% hold on

% figure(7)

% plot(M,Nux,Cl(i),LineWidth=2);

% xlabel('M');

% ylabel("$Re_{x}^{-(1/2)}Nu_{x}$",Interpreter="latex");

% hold on

% figure(8)

% plot(M,Shx,Cl(i),LineWidth=2);

% xlabel('M');

% ylabel("$Re_{x}^{-(1/2)}Sh_{x}$",Interpreter="latex");

% hold on

figure(1);