Getting error while trying to solve boundary condition equations.

I'm currently working on solving a boundary value problem (BVP) represented by a system of ordinary differential equations in MATLAB. However, I've encountered some errors during the process, and I'm seeking assistance in resolving them. Below is the code snippet

close all; clear all; clc;

K1=0; eta1=0.2; eta2=0.7; lambda1=0.5; Fr=0.5; Rd=0.5; Sc=0.8; Pr=1; omega=0.5; Nr=0.5; B1=0.5; Nt=0.4; beta1=0.3; m1=0.4; Pe=0.3; E2=0.2; Lb=0.4; Rb=0.6; alpha=0.3; r1=0.3; Kp=0.3;
w=1; Nb=0.5; B2=0.5;
% fprintf('Working\n');
for M1 = 0:0.1:1
    V1(w)= M1;
   
    j=1;
% fprintf('Working\n');
% global khp;
% khp =@ (x,y) [y(2);
%     y(3);
%     (-y(1)*y(3)+(1+Fr)*(y(2)^2)+Kp*y(2)+M1*y(2)-lambda1*(y(5)-Nr*y(8)-Rb*y(11)))/(1+K1-eta1*K1*(y(3)^2));
%     y(6);
%   -(Pr*y(6)*y(1)+Pr*Nb*y(6)*y(9)+Pr*Nt*((y(6))^2)+alpha*Pr*y(5))/(1+Rd*4/3);
%    y(9);
%     (Sc*beta1*(((1+y(5)*r1))^M1)*y(8)*exp(-E2/(1+r1*y(5)))-Sc*y(1)*y(9)-(Nt/Nb)*y(7));
%     y(12);
%     (Pe*(y(11)+omega)*y(10)-y(12)*(Lb*y(1)-Pe*y(9)))];

 global khp;
khp =@ (x,y) [y(2);
    y(3);
    ((-y(1)*y(3))+((1+Fr)*(y(2)^2))+(Kp*y(2))+(M1*y(2))-(lambda1*(y(4)-Nr*y(6)-Rb*y(8))))/((1+K1)-(eta1*K1*(y(3)^2)));
    y(5);
  (-((Pr*y(5)*y(1))+(Pr*Nb*y(5)*y(7))+(Pr*Nt*((y(5))^2))+(alpha*Pr*y(4)))/(1+(Rd*4/3)));
    y(7);
    (Sc*beta1*(((1+(y(4)*r1)))^m1)*y(6)*exp(-E2/(1+(r1*y(4))))-(Sc*y(1)*y(7))-((Nt/Nb)*(-((Pr*y(5)*y(1))+(Pr*Nb*y(5)*y(7))+(Pr*Nt*((y(5))^2))+(alpha*Pr*y(4)))/(1+(Rd*4/3)))));
   y(9);
  (Pe*(y(8)+omega)*((Sc*beta1*(((1+(y(4)*r1)))^m1)*y(6)*exp(-E2/(1+(r1*y(4))))-(Sc*y(1)*y(7))-((Nt/Nb)*(-((Pr*y(5)*y(1))+(Pr*Nb*y(5)*y(7))+(Pr*Nt*((y(5))^2))+(alpha*Pr*y(4)))/(1+(Rd*4/3))))))+((y(9))*((Pe*y(6))-(Lb*y(1)))))];


% fprintf('Working\n');
init=bvpinit(linspace(0,9,100),[1/2, 0, 0, 1, 0, 0, 1, 0, 0]);
% fprintf('Working\n');
bc_output = bc_bvp(sol1);
sol1=bvp4c(khp,bc_output,init);
% fprintf('Working\n');
x=linspace(0,9,100);
y1=deval(sol1,x);

Cf(w,j)= ((1+K1)*y1(3,1)-((K1*eta1)/3)*y1(3,1)^(3));
Nu(w,j)=-(1+Rd*4/3)*y1(5,1);
Sh(w,j)=-(y1(7,1));
Nn(w,j)=-(y1(9,1));

w=w+1;
end
%
M1 = 0:0.1:1;
plot(M1,Cf)
function res = bc_bvp(yl,yr) % boundary conditions
eta1=0.2; eta2=0.7;K1=0;B1=0.5; B2=0.5;
%global khp;
%init=bvpinit(linspace(0,9,100),[1/2, 0, 0, 1, 0, 0, 1, 0, 0, 0,0,0]);
bc_output = bc_bvp(sol1);
sol1=bvp4c(khp,bc_output,init);
%y1=deval(sol1,x);
% fprintf("%f\n",length(yr));
%fprintf("%f\n",yr(11));
res=[yl(1); yl(2)-1-(eta2)*((1+K1)*yl(3)-(eta1)*K1*(yl(3)^3)); yl(5)+B1*(1-yl(4)); yl(7)+B2*(1-yl(6)); yl(9)-1; yr(2); yr(4); yr(6); yr(9);];
% fprintf("%f\n",length(res));
end

looking for some guidance or suggestions .

Answer 1

[Pavan Sahith]

It appears that you are encountering an issue while trying to solve a boundary value problem (BVP) using MATLAB. The error you're experiencing is due to the use of the variable ‘sol1’ before it has been properly defined. To address this problem, you need to ensure that ‘sol1’ is defined by calling “bvp4c” before attempting to access or use it.

Issue might be solved by passing “bc_bvp” as a function handle to “bvp4c” instead of using the 'sol1' before defining.

In MATLAB, the “bvp4c” function is designed to solve BVPs and requires three key arguments:

1. The function that defines the differential equations (‘khp’ in your case).
2. The function that specifies the boundary conditions (‘bc_bvp’ in your case).
3. An initial guess for the solution structure (‘init’ in your case).

You can pass the boundary condition function “bc_bvp” as a function handle using the @ symbol. The correct way to call ‘bvp4c” and define ‘sol1’ would be as follows:

sol1 = bvp4c(khp, @bc_bvp, init); 

For additional guidance and examples on how to use “bvp4c”, including how to pass parameters to the function, you can refer to the MathWorks documentation.

https://www.mathworks.com/help/matlab/ref/bvp4c.html

Hope this will help you in solving the error.

Answer 2

[Torsten]

For M1 = 0-0.4, bvp4c cannot find a solution. So I started with M1 = 0.5.

close all; clear all; clc;

K1=0; eta1=0.2; eta2=0.7; lambda1=0.5; Fr=0.5; Rd=0.5; Sc=0.8; Pr=1; omega=0.5; Nr=0.5; B1=0.5; Nt=0.4; beta1=0.3; m1=0.4; Pe=0.3; E2=0.2; Lb=0.4; Rb=0.6; alpha=0.3; r1=0.3; Kp=0.3;
w=1; Nb=0.5; B2=0.5;
% fprintf('Working\n');
for M1 = 0.5:0.1:1
    V1(w)= M1;
   
    %j=1;
% fprintf('Working\n');
% global khp;
% khp =@ (x,y) [y(2);
%     y(3);
%     (-y(1)*y(3)+(1+Fr)*(y(2)^2)+Kp*y(2)+M1*y(2)-lambda1*(y(5)-Nr*y(8)-Rb*y(11)))/(1+K1-eta1*K1*(y(3)^2));
%     y(6);
%   -(Pr*y(6)*y(1)+Pr*Nb*y(6)*y(9)+Pr*Nt*((y(6))^2)+alpha*Pr*y(5))/(1+Rd*4/3);
%    y(9);
%     (Sc*beta1*(((1+y(5)*r1))^M1)*y(8)*exp(-E2/(1+r1*y(5)))-Sc*y(1)*y(9)-(Nt/Nb)*y(7));
%     y(12);
%     (Pe*(y(11)+omega)*y(10)-y(12)*(Lb*y(1)-Pe*y(9)))];

 %global khp;
khp =@ (x,y) [y(2);
    y(3);
    ((-y(1)*y(3))+((1+Fr)*(y(2)^2))+(Kp*y(2))+(M1*y(2))-(lambda1*(y(4)-Nr*y(6)-Rb*y(8))))/((1+K1)-(eta1*K1*(y(3)^2)));
    y(5);
  (-((Pr*y(5)*y(1))+(Pr*Nb*y(5)*y(7))+(Pr*Nt*((y(5))^2))+(alpha*Pr*y(4)))/(1+(Rd*4/3)));
    y(7);
    (Sc*beta1*(((1+(y(4)*r1)))^m1)*y(6)*exp(-E2/(1+(r1*y(4))))-(Sc*y(1)*y(7))-((Nt/Nb)*(-((Pr*y(5)*y(1))+(Pr*Nb*y(5)*y(7))+(Pr*Nt*((y(5))^2))+(alpha*Pr*y(4)))/(1+(Rd*4/3)))));
   y(9);
  (Pe*(y(8)+omega)*((Sc*beta1*(((1+(y(4)*r1)))^m1)*y(6)*exp(-E2/(1+(r1*y(4))))-(Sc*y(1)*y(7))-((Nt/Nb)*(-((Pr*y(5)*y(1))+(Pr*Nb*y(5)*y(7))+(Pr*Nt*((y(5))^2))+(alpha*Pr*y(4)))/(1+(Rd*4/3))))))+((y(9))*((Pe*y(6))-(Lb*y(1)))))];


% fprintf('Working\n');
init=bvpinit(linspace(0,9,100),[1/2, 0, 0, 1, 0, 0, 1, 0, 0]);
% fprintf('Working\n');
%bc_output = bc_bvp(sol1);
sol1=bvp4c(khp,@bc_bvp,init);

% fprintf('Working\n');
x=linspace(0,9,100);
y1=deval(sol1,x);

Cf(w)= ((1+K1)*y1(3,1)-((K1*eta1)/3)*y1(3,1)^(3));
Nu(w)=-(1+Rd*4/3)*y1(5,1);
Sh(w)=-(y1(7,1));
Nn(w)=-(y1(9,1));

w=w+1;
end
%
M1 = 0.5:0.1:1;
plot(M1,Cf)

function res = bc_bvp(yl,yr) % boundary conditions
eta1=0.2; eta2=0.7;K1=0;B1=0.5; B2=0.5;
%global khp;
%init=bvpinit(linspace(0,9,100),[1/2, 0, 0, 1, 0, 0, 1, 0, 0, 0,0,0]);
%bc_output = bc_bvp(sol1);
%sol1=bvp4c(khp,bc_output,init);
%y1=deval(sol1,x);
% fprintf("%f\n",length(yr));
%fprintf("%f\n",yr(11));
res=[yl(1) 
    yl(2)-1-(eta2)*((1+K1)*yl(3)-(eta1)*K1*(yl(3)^3)) 
    yl(5)+B1*(1-yl(4)) 
    yl(7)+B2*(1-yl(6)) 
    yl(9)-1 
    yr(2) 
    yr(4) 
    yr(6) 
    yr(9)];
% fprintf("%f\n",length(res));
end