How to display Skinfriction and Nusslet output in marix form and table in workspace, error in line 40

 %bvp4c_sys5.m       5 coupled 1st order ODEs  in a 2-point BVP
clc;clear;close;
 A1=2.5;
A2=6.5;
m=3.0;
beta=2.2569438;
S=0.2;
K = 0.5;

eta =0;
gamma = 0.1;
Pr=7.56;
kf = 100;
ks = 0.1;
g1 = 0.001;
g2 = 0.001;
phi=0.02;
rf=0.00001;
rs=0.01;
Q=2.0;
ss=0.000001;
sf=0.1;
ps=0.000001;
pf=0.1;
Ec=1.0;
knf = kf*(((((ks+(m-1)*kf))+((m-1)*(ks-kf)*phi)))/(ks+((m-1)*kf)-(ks-kf)*phi));
e1 = (1+A1*phi+A2*((phi)^2))/(1-phi+phi*(rs/rf));
e2 = (knf/kf)/(1-phi+phi*(ps/pf));
e3 = (1-phi+phi*(ss/sf))/(1-phi+phi*(rs/rf));

 x=linspace(0,1,5);
 % initial mesh.
 yinit=[g1 g2 (-1/K)*g2 0 -((kf/knf)*gamma*(1-g2))];
 % Guess  vector for the solution
 M=linspace(0,1,100);
for S=0.2:0.2:0.6
 ode = @(x,y) [ y(2); y(3); (((e3/e1)*M*y(2))-((1/e1)*((y(1)*y(3))-y(2)^2-(S*(y(2)+(eta/2)*y(3)))))); y(5); ((-(e1/e2)*Pr*Ec*y(3)^2)-((Pr/e2)*(((y(1)*y(5))-2*y(4)*y(2)-(Pr/e2)*Q*y(4))+(Pr/e2)*((S/2)*(3*y(4)+(eta*y(5))))))) ];
 bc= @(ya,yb) [ ya(2)-1-K*ya(3); ya(5)+((kf/knf)*gamma*(1-ya(4))); yb(1)-((S*beta)/2); yb(3); yb(5) ]; % Boundary conditions
 solinit = bvpinit(x,yinit) % initial guess 
  sol=bvp4c(ode, bc, solinit); 
 y = deval(sol,x);% Evaluation of y at x
 [M' y'];
 sk=(1.0+A1*phi+A2*((phi)^2)*y(3, 1));
 %X=M'; Y=y(1,:)'; dYdX=y(2,:)'; sk=((1.0+A1*phi+A2*((phi)^2)*y(3, 1)))';  
 %A = table(X,Y,dYdX,(1.0+A1*phi+A2*((phi)^2)*y(3, 1)))
  figure;
plot(M,y(5,1),'b','LineWidth',1.2)
ylabel('{\theta}({\eta})')
xlabel('{\eta}')   
 plot(M, -(1.0+A1*phi+A2*((phi)^2)*y(3)),'linewidth',2); 
 %xlim=[0 1];
 %ylim=[-2 1];
 hold on
 fprintf('Cfx = %7.9f\n',-y(3,1)*(1.0+A1*phi+A2*((phi)^2)));
   
end
 figure;
%plot(x,y(5,1),'b','LineWidth',1.2)
%ylabel('{\theta}({\eta})')
%xlabel('{\eta}')   
grid; legend('S=0.2','S=0.4','S=0.6');
 xlabel('\eta'); ylabel('f_1(\eta)');

Answer 1

[Shivam]

Hi Lakshana,

I understand that you get an error while trying to find skin friction while executing the provided code.

Upon investigating the error, I see that the M vector, which has a dimension of 1x100, is used incorrectly inside the differential equation definition (line 37). Here, the expression (e3/e1)*M*y(2) generates an array with dimensions 1x100 while other row elements are scalars. This dimension inconsistency while attempting to concatenate the rows leads to the error. You can note that to concatenate the arrays or matrices vertically, they must have the same number of columns. Additionally, there is a misplacement in the sequence of plot commands.

You can refer to the modified code, which is adjusted to correctly loop over different values of S and plot the results for each case in a single figure. Also, M is not directly accessed within the differential function ode.

% .
% Parameters Initialization
% ..
x = linspace(0, 1, 5); % Initial mesh.
yinit = [g1, g2, (-1/K)*g2, 0, -((kf/knf)*gamma*(1-g2))]; % Guess vector for the solution.

figure;
colors = ['r', 'g', 'b']; % Color for 3 S value plot.
cIndex = 1;

for S = 0.2:0.2:0.6
    ode = @(x,y) [y(2); y(3); (((e3/e1)*y(2)) - ((1/e1)*((y(1)*y(3)) - y(2)^2 - (S*(y(2) + (eta/2)*y(3)))))); y(5); ((-(e1/e2)*Pr*Ec*y(3)^2) - ((Pr/e2)*(((y(1)*y(5)) - 2*y(4)*y(2) - (Pr/e2)*Q*y(4)) + (Pr/e2)*((S/2)*(3*y(4) + (eta*y(5)))))))];
    bc = @(ya,yb) [ya(2)-1-K*ya(3); ya(5)+((kf/knf)*gamma*(1-ya(4))); yb(1)-((S*beta)/2); yb(3); yb(5)]; % Boundary conditions.
    solinit = bvpinit(x, yinit); % Initial guess.
    sol = bvp4c(ode, bc, solinit);
    y = deval(sol, x); % Evaluation of y at x.

    % Plotting
    plot(x, y(5,:), colors(cIndex), 'LineWidth', 1.2); hold on;
    cIndex = cIndex + 1; % Move to next color.

    fprintf('Cfx = %7.9f for S = %4.2f\n', -y(3,1)*(1.0 + A1*phi + A2*(phi^2)), S);
end

Cfx = 1.528551677 for S = 0.20
Cfx = 1.751184435 for S = 0.40
Cfx = 1.848326934 for S = 0.60

xlabel('\eta'); ylabel('f_1(\eta)');
legend('S=0.2', 'S=0.4', 'S=0.6'); 
grid on;
title('Plot of f_1(\eta) for different values of S');
hold off;

I trust it addresses your issue.

Thanks