In my code, I am getting plot on lambda=1, but rest plots are not reflecting. Please help to get all figures.
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
Accepted Answer
The results for the last 3 values of lambda are NaN (see above) - thus they are not plotted.
2 Comments
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
Reshu
on 15 Feb 2026
More Answers (0)
Categories
Find more on Fortran with MATLAB in Help Center and File Exchange
