Looking at the screenshots, there are multiple values of M used to plot the graphs for F but your function only uses a single value (M = 0.2). To get multiple lines in the graph, you need to modify the function to accept different values of M, or accept all the values of M as a vector and use a for-loop to plot the result for each value.
Here is an example where all the values of M are passed as a vector and then plotted in the same graph:
conedisk_multiple_M([0.2 0.3 0.4])
function conedisk_multiple_M(M_values)
% Parameter values
A1 = 1.10629;
A2 = 1.15;
A3 = 1.2;
A4 = 1.1;
A5 = 1.1;
A6 = 1.1;
Grt = 5;
Pr = 0.71;
R = 0.2;
Ec = 0.1;
Q = 0.1;
Rew = 12;
Red = -12;
n = -1;
g1 = 0;
g2 = 0;
g3 = 0;
g4 = 0;
inf = 1;
% Set solver options to increase the maximum number of mesh points
options = bvpset('RelTol', 1e-5, 'AbsTol', 1e-7, 'NMax', 5000);
% Loop over each value of M
for M = M_values
% Defining parameters
solinit = bvpinit(linspace(0, 1, 200), [0 g1 g2 0 Rew g3 1 g4]);
sol = bvp4c(@bvp2D, @bc2D, solinit);
x = sol.x;
y = sol.y;
% Plotting of the velocity
figure(1)
plot(x, y(1,:), 'DisplayName', sprintf('M = %.2f', M), 'linewidth', 1)
hold on
xlabel('\eta', 'fontweight', 'bold', 'fontsize', 16)
ylabel('F(\eta)', 'fontweight', 'bold', 'fontsize', 16)
% Similarly plot G and theta
end
% Add legends to the plots
figure(1)
legend show
hold off
% Residual of the boundary conditions
function res = bc2D(y0, yinf)
res = [y0(1);yinf(1); y0(4); yinf(4);y0(5)-Rew; yinf(5)-Red;y0(7)-1; yinf(7)];
end
% System of First Order ODEs
function dydx = bvp2D(eta, y)
yy1 = (-(1+eta^2)*((10*A1*eta+A2*eta*y(1)-A2*y(4))*y(3))-3*(2*A1+7*A1*(eta^2)+A2*(2*eta^2+1)*y(1)-A2*eta*y(4))*y(2)-2*A2*y(5)*y(6)-3*(A1+A2*y(1))*y(4)+A3*M*y(2)-A4*Grt*y(8))/(A1*(1+eta^2)^2);
yy2 = eta*y(2);
yy3 = (-3*A1*eta*y(6)-A2*(eta*y(1)-y(4))*y(6)+A3*M*y(5))/(A1*(eta^2+1));
yy4 = ((A5+4/3*R)*((eta*((2*n)-1))*y(8)+n^2*y(7))-A3*Pr*M*Ec*(y(1)^2+y(5)^2)-Pr*Q*y(7)+A6*Pr*(y(4)*y(8)+n*y(1)*y(7)-eta*y(1)*y(8)))/(A5+4/3*R)*(1+eta^2);
dydx = [y(2); y(3); yy1; yy2; y(6); yy3; y(8); yy4];
end
end
The current results do not match the graph shown in the screenshots, suggesting that there are some issues with the equations you are using to solve the problem. I would suggest checking the equations again to make sure they are correct. Once you have fixed the equations, you should get the multiple graphs as expected.
If you also want to plot the second F graph, you need to modify your function to accept different values for Rew and Red since the second graph uses 
, 
instead of 
, 
. You can then use the subplot function to create a tiled plot. For example:
, instead of , . You can then use the subplot function to create a tiled plot. For example:conedisk_subplot([0.2 0.3 0.4], [12 0], [0 12])
function conedisk_subplot(M_values, Rew_values, Red_values)
% Parameter values
A1 = 1.10629;
A2 = 1.15;
A3 = 1.2;
A4 = 1.1;
A5 = 1.1;
A6 = 1.1;
Grt = 5;
Pr = 0.71;
R = 0.2;
Ec = 0.1;
Q = 0.1;
n = -1;
g1 = 0;
g2 = 0;
g3 = 0;
g4 = 0;
inf = 1;
% Set solver options to increase the maximum number of mesh points
options = bvpset('RelTol', 1e-5, 'AbsTol', 1e-7, 'NMax', 5000);
% Determine the number of subplots needed
num_plots = length(Rew_values);
num_cols = ceil(sqrt(num_plots));
num_rows = ceil(num_plots / num_cols);
% Create a new figure for the plots
figure('Name', 'Velocity Plots for Different M Values')
% Loop over each pair of Rew and Red
for i = 1:num_plots
Rew = Rew_values(i);
Red = Red_values(i);
% Create a subplot for each pair of Rew and Red
subplot(num_rows, num_cols, i)
hold on
% Loop over each value of M
for M = M_values
% Defining parameters
solinit = bvpinit(linspace(0, 1, 200), [0 g1 g2 0 Rew g3 1 g4]);
sol = bvp4c(@bvp2D, @bc2D, solinit);
x = sol.x;
y = sol.y;
% Plotting of the velocity
plot(x, y(1,:), 'DisplayName', sprintf('M = %.2f', M), 'linewidth', 1)
end
% Customize subplot
title(sprintf('$Re_\\omega = %.0f, Re_\\Omega = %.0f$', Rew, Red), 'Interpreter', 'latex')
xlabel('\eta', 'fontweight', 'bold', 'fontsize', 12)
ylabel('F(\eta)', 'fontweight', 'bold', 'fontsize', 12)
legend show
grid on
hold off
end
% Residual of the boundary conditions
function res = bc2D(y0, yinf)
res = [y0(1);yinf(1); y0(4); yinf(4);y0(5)-Rew; yinf(5)-Red;y0(7)-1; yinf(7)];
end
% System of First Order ODEs
function dydx = bvp2D(eta, y)
yy1 = (-(1+eta^2)*((10*A1*eta+A2*eta*y(1)-A2*y(4))*y(3))-3*(2*A1+7*A1*(eta^2)+A2*(2*eta^2+1)*y(1)-A2*eta*y(4))*y(2)-2*A2*y(5)*y(6)-3*(A1+A2*y(1))*y(4)+A3*M*y(2)-A4*Grt*y(8))/(A1*(1+eta^2)^2);
yy2 = eta*y(2);
yy3 = (-3*A1*eta*y(6)-A2*(eta*y(1)-y(4))*y(6)+A3*M*y(5))/(A1*(eta^2+1));
yy4 = ((A5+4/3*R)*((eta*((2*n)-1))*y(8)+n^2*y(7))-A3*Pr*M*Ec*(y(1)^2+y(5)^2)-Pr*Q*y(7)+A6*Pr*(y(4)*y(8)+n*y(1)*y(7)-eta*y(1)*y(8)))/(A5+4/3*R)*(1+eta^2);
dydx = [y(2); y(3); yy1; yy2; y(6); yy3; y(8); yy4];
end
end
Similarly, the graph of G is plotted using different values of 
and you can use a similar approach as above.
and you can use a similar approach as above.Please refer to the following resources for more information:
- subplot documentation - https://www.mathworks.com/help/matlab/ref/subplot.html
- hold documentation - https://www.mathworks.com/help/matlab/ref/hold.html
- for-loop documentation - https://www.mathworks.com/help/matlab/ref/for.html
Hope this helps!
