Try this (change the ii to satisfy whatever range of N you would like to observe):
for ii = 0:0.1:1
Entropy(ii)
end
Where I modified your function to take an input:
function [] = Entropy(M)
clc
%M=0.0;
Pr=1.0; %Prandtl
Ec=0.5;
Bi=0.5;
Omega=1.0;
Re=5.0;
Br=1.0;
% Initial values
sol = bvpinit(linspace(0,8,30),[1 0 0 0 0 ]);
% solution in structure form
sol1 = bvp4c(@bvpexam2, @bcexam2, sol);
% x values
x1 = sol1.x;
% y values in row form (y, y', y'', theta, theta');
y1 = sol1.y;
figure (1)
plot(x1, y1(4, :),'-R','linewidth',1,'Linesmoothing','on');
hold on
grid on
% Here I define residual of boundary conditions
function res = bcexam2(y0, yinf)
res = [y0(1); y0(2)-1; y0(5)+Bi*(1-y0(4)); yinf(2); yinf(4)];
end
% First order ODEs are define here
function ysol = bvpexam2(x,y)
yy1=-2*y(1)*y(3)+y(2)*y(2)+M*y(2);
yy2=-2*Pr*y(1)*y(5)-Pr*Ec*y(3)*y(3)-M*Pr*Ec*y(2)*y(2);
Ns=Re*y(5)*y(5)+Re*(Br/Omega)*y(3)*y(3)+Re*(Br/Omega)*M*y(2)*y(2);
Be=Re*y(5)*y(5)/(Re*y(5)*y(5)+Re*(Br/Omega)*y(3)*y(3)+Re*(Br/Omega)*M*y(2)*y(2));
ysol = [y(2);y(3);yy1; y(5);yy2];
end
end
