In my code, I am getting plot on lambda=1, but rest plots are not reflecting. Please help to get all figures.
35 views (last 30 days)
Show older comments
clc; clear; close all;
%% Parameters (from paper)
Pr = 0.71;
%%lambda = 0.2;
N = 0.5;
zeta=0.1;
cosdlt = 0.3;
epsi = 1;
theta_w = 1.5;
S = 0.5;
M = 0;
%% Zeta values
lambda_list = [1 3 5 7];
%% Numerical settings
z0 = 0;
zinf = 5;
h = 0.01;
figure(1); hold on; % f'
figure(2); hold on; % theta
% Initial guess for zeta = 0 (important)
s_prev = [0.75; -0.6]; % [f'(0), theta'(0)]
for m = 1:length(lambda_list)
lambda = lambda_list(m);
par = [Pr lambda N zeta cosdlt epsi theta_w S];
% Continuation-based guesses
s1 = s_prev;
s2 = s_prev + [0.05; -0.05];
% Secant shooting
for k = 1:30
F1 = shoot(s1,par,z0,zinf,h);
F2 = shoot(s2,par,z0,zinf,h);
s = s2 - (s2-s1).*F2./(F2-F1);
if norm(s-s2) < 1e-6
break
end
s1 = s2;
s2 = s;
end
sol = s;
s_prev = sol; % continuation step
% Final integration
y0 = [0; sol(1); -2*S; 1; sol(2)];
[z,y] = rkf45(@(z,y) odesys(z,y,par), [z0 zinf], y0, h);
y
% ---- Plots ----
figure(1)
plot(z,y(:,2),'LineWidth',2)
figure(2)
plot(z,y(:,4),'LineWidth',2)
end
%% Formatting
figure(1)
xlabel('\eta'); ylabel('f''(\eta)')
legend('\lambda=1','\lambda=3','\lambda=5','\lambda=7','Location','northeast')
grid on; box on;
figure(2)
xlabel('\eta'); ylabel('\theta(\eta)')
legend('\lambda=1','\lambda=3','\lambda=5','\lambda=7','Location','northeast')
grid on; box on;
% Save solution for streamlines
eta = z;
f = y(:,1);
save flow_solution.mat eta f
function dydz = odesys(z,y,par)
Pr = par(1);
lambda = par(2);
N = par(3);
zeta = par(4);
cosdlt = par(5);
epsi = par(6);
theta_w = par(7);
M = par(8);
f = y(1);
fp = y(2);
fpp = y(3);
th = y(4);
thp = y(5);
% Governing equations
fppp = fp^2 - 2*f*fpp + lambda*fp - zeta*th*cosdlt -M^2*fp;
thpp = ( ...
-2*Pr*f*thp + 2*Pr*fp*th ...
+ Pr*epsi*(th + 1/(theta_w-1))*fp ...
) / (1+N);
dydz = [
fp
fpp
fppp
thp
thpp
];
end
function [t,y] = rkf45(odefun,tspan,y0,h)
t = (tspan(1):h:tspan(2))';
n = length(t);
y = zeros(n,length(y0));
y(1,:) = y0';
for i = 1:n-1
ti = t(i);
yi = y(i,:)';
k1 = h*odefun(ti,yi);
k2 = h*odefun(ti+h/4, yi+k1/4);
k3 = h*odefun(ti+3*h/8, yi+3*k1/32+9*k2/32);
k4 = h*odefun(ti+12*h/13, yi+1932*k1/2197 ...
-7200*k2/2197 +7296*k3/2197);
k5 = h*odefun(ti+h, yi+439*k1/216 ...
-8*k2 +3680*k3/513 -845*k4/4104);
k6 = h*odefun(ti+h/2, yi ...
-8*k1/27 +2*k2 -3544*k3/2565 ...
+1859*k4/4104 -11*k5/40);
y(i+1,:) = (yi + ...
16*k1/135 +6656*k3/12825 ...
+28561*k4/56430 -9*k5/50 ...
+2*k6/55)';
end
end
function F = shoot(s,par,z0,zinf,h)
S = par(end);
% Initial conditions
y0 = [ ...
0; % f(0)
s(1); % f'(0) (guessed)
-2*S; % f''(0) (given)
1; % theta(0)
s(2) % theta'(0) (guessed)
];
[z,y] = rkf45(@(z,y) odesys(z,y,par), [z0 zinf], y0, h);
% Boundary residuals
F = [y(end,2); y(end,4)];
% Safety against divergence
if any(isnan(F)) || any(isinf(F))
F = [1;1];
end
end
2 Comments
Accepted Answer
Torsten
on 14 Feb 2026 at 12:34
Edited: Torsten
on 14 Feb 2026 at 12:36
The results for the last 3 values of lambda are NaN (see above) - thus they are not plotted.
2 Comments
Torsten
on 14 Feb 2026 at 15:22
Edited: Torsten
on 15 Feb 2026 at 17:18
I noticed that you solve a boundary value problem. Why don't you use MATLAB's "bvp4c" ?
clc; clear; close all;
%% Parameters (from paper)
Pr = 0.71;
%%lambda = 0.2;
N = 0.5;
zeta=0.1;
cosdlt = 0.3;
epsi = 1;
theta_w = 1.5;
S = 0.5;
M = 0;
%% Zeta values
lambda_list = [1 3 5 7];
%% Numerical settings
z0 = 0;
zinf = 5;
xmesh = linspace(z0,zinf,100);
solinit = bvpinit(xmesh, [0 0 0 0 0]);
figure(1)
hax1=axes;
hold on
figure(2)
hax2=axes;
hold on
for i = 1:numel(lambda_list)
lambda = lambda_list(i);
par = [Pr lambda N zeta cosdlt epsi theta_w S];
sol = bvp4c(@(x,y)odesys(x,y,par), @(ya,yb)odebc(ya,yb,par), solinit);
plot(hax1,sol.x,sol.y(2,:),'LineWidth',2)
plot(hax2,sol.x,sol.y(4,:),'LineWidth',2)
end
figure(1)
xlabel('\eta'); ylabel('f''(\eta)')
legend('\lambda=1','\lambda=3','\lambda=5','\lambda=7','Location','northeast')
grid on; box on;
hold off
figure(2)
xlabel('\eta'); ylabel('\theta(\eta)')
legend('\lambda=1','\lambda=3','\lambda=5','\lambda=7','Location','northeast')
grid on; box on;
hold off
%--------------------------------
function res = odebc(ya,yb,par) % boundary conditions
S = par(end);
res = [ya(1)
yb(2)
ya(3)+2.0*S
ya(4)-1.0
yb(4)
];
end
%--------------------------------
function dydz = odesys(z,y,par)
Pr = par(1);
lambda = par(2);
N = par(3);
zeta = par(4);
cosdlt = par(5);
epsi = par(6);
theta_w = par(7);
M = par(8);
f = y(1);
fp = y(2);
fpp = y(3);
th = y(4);
thp = y(5);
% Governing equations
fppp = fp^2 - 2*f*fpp + lambda*fp - zeta*th*cosdlt -M^2*fp;
thpp = ( ...
-2*Pr*f*thp + 2*Pr*fp*th ...
+ Pr*epsi*(th + 1/(theta_w-1))*fp ...
) / (1+N);
dydz = [
fp
fpp
fppp
thp
thpp
];
end
More Answers (0)
See Also
Categories
Find more on Boundary Value Problems in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
You can also select a web site from the following list
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)