sol = bvp4c (OdeBVP, OdeBC, solinit, options);

ne the boundary conditions

function res = OdeBc (ya, yb, A, s, B, lambda)

global A s B lambda

res= [ya(1)-s;

ya(2)-lambda-A*ya(3);

ya(4)-1-B*ya(5);

yb(2);

yb(4)];

end

% setting the initial guess for first solution

function v = OdeInit1(x,A,s,lambda)

global A s lambda

v=[s+0.56

0

0

0

0];

end

% setting the initial guess for second solution

function v1 =OdeInit2(x, A, s)

global A s

v1 = [exp(-x)

exp(-x)

-exp(-x)

-exp(-x)

-exp(-x)];

end

end

[Torsten]

solinit = bvpinit (linspace (etaMin, etaMax1, stepsize1),@(x)OdeInit1(x,A,s,lambda));
sol = bvp4c (@(x,y)OdeBVP(x,y,Pr,D1,Kh,H4,H3,beta), @(ya,yb)OdeBC(ya, yb, A, s, B, lambda), solinit, options);

Same for your second call to bvp4c.

And define the parameters before calling bvp4c.

And remove the global statements inside the functions.

[KIRAN Sajjanshetty]

can u change directly in the code and send it here

i have tried for all possible ways but it is not running showing error

[Torsten]

Include the modified code (with all parameters you use) and the error message here.

[KIRAN Sajjanshetty]

Unrecognized function or variable 's'.

Error in slipflow>@(x)OdeInit1(x,A,s,lambda) (line 59)

solinit = bvpinit (linspace (etaMin, etaMax1, stepsize1),@(x)OdeInit1(x,A,s,lambda));

Error in bvpinit (line 104)

w = feval(v,x(1),extraArgs{:});

Error in slipflow (line 59)

solinit = bvpinit (linspace (etaMin, etaMax1, stepsize1),@(x)OdeInit1(x,A,s,lambda));

[Torsten]

Seems you didn't define "s" and gave it a value before calling bvpinit.

[Waseef]

sir this my code for generating Nusselt number plot but i got an empty plot can you correct it for me. I want to plot it aganist Phi or aa

function waslipflow

format long g

%Define all parameters

% Boundary layer thickness & stepsize

etaMin = 0.01; etaMax1 = 10; etaMax2 = 15; stepsize1 = etaMax1;

stepsize2 = etaMax2;

% Input for the parameters

Phi=0.01;%input ('Input the value of phi = ');

R =0; %input ('Input the value of Radiation = ');

aa = 10;%input ('Input the value of aa = ');

M= 0.0;%input('magentic parameter M =');

pm=0.0;%input('porosity =');

Q=0.0; %input('Heat sourse parameter');

sigmaf=0.05;

sigmas=59600000;

n=3;%input ('Input the value of n = ');

ks= 400; kf= 0.613; Rhos= 8933; Rhof= 997.1;

Cps= 385; Cpf= 4179; tt=Rhos*Cps; kk=Rhof*Cpf;

Pr = 10; % input ('Input the Prandtl number = ');

Ec =0.0;%input ('Input Eckret for velociy exponent parameter = ');

spn=(3*Phi*((sigmas/sigmaf)-1))/(((sigmas/sigmaf)+2)-((sigmas/sigmaf)-1)*Phi);

sigmanf=sigmaf*(1+spn);

phi1=(1-Phi)^2.5;

phi2=1-Phi+Phi*(Rhos/Rhof);

A=phi1*phi2;

B=phi1*phi2*(2*aa/(aa+1));

C=(2/(aa+1))*((phi1*M*(sigmanf/sigmaf))+pm);

total=sigmas/sigmaf;

phi3=1+(3*(total-1)*Phi/((total+2)-(total-1)*Phi));

phi4=(ks+(n-1)*kf+(n-1)*(ks-kf)*Phi)/(ks+(n-1)*kf+(kf-ks)*Phi);

phi5=1-Phi+Phi*(Cps/Cpf);

phi5=1-Phi+Phi*(tt/kk);

A1=phi4+R;

A2=Pr*phi5*(4*aa/(aa+1));

A3=Pr*phi5;

A4=(2/(aa+1))*Pr*Q;

A5=Ec*Pr/phi1;

A6=(2/(aa+1))*Pr*Ec*M*(sigmanf/sigmaf);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%First solution %%%%%%%%%%%%%%%%%%%%%%

options = bvpset('stats','off','RelTol',1e-6);

solinit = bvpinit (linspace (etaMin, etaMax1, stepsize1),@(x)OdeInit1(x,A));

sol = bvp4c (@(x,y)OdeBVP(x,y,A,B,C,A1,A2,A3,A4,A5,A6), @(ya,yb)OdeBc(ya, yb), solinit, options);

eta = linspace (etaMin, etaMax1, stepsize1);

y= deval (sol,eta);

% Calculate the Nusselt number

%Nu = -phi4*sqrt((aa+1/2)) * sol.y(5,1);

% Plot the Nusselt number

figure(3)

plot(eta, -phi4*sqrt((aa+1/2)) * sol.y(5,1), '-')

xlabel('\eta')

ylabel('Nu')

xlim([0 8])

ylim([0 1])

title('Nusselt Number Profile')

% figure(1) %velocity profile

% plot(sol.x,sol.y(2,:),'-')

% xlabel('\eta')

% ylabel('f`(\eta)')

% hold on

% figure(2) %temperature profile

% plot(sol.x,sol.y(4,:),'-')

% xlabel('\eta')

% ylabel('\theta(\eta)')

% xlim([0 8])

% ylim([0 1])

% hold on

% saving the out put in text file for first solution

descris =[sol.x; sol.y];

%save 'sliphybrid_upper.txt' descris -ascii

% Displaying the output for first solution

% fprintf('\n First solution:\n');

% fprintf('f"(0)=%7.9f\n',y(3,:)); % reduced skin friction

% fprintf('-theta(0)=%7.9f\n',-y(5,:)); %reduced local nusselt number

% fprintf('Cfx=%7.9f\n',B*(y(3,:))); % skin friction

% fprintf('Nux=%7.9f\n',-A*y(5,:)); % local nusselt number

% fprintf('\n');

end

% Define the ODE function

function f = OdeBVP(~,y,A,B,C,A1,A2,A3,A4,A5,A6)

f =[y(2); y(3);B*y(2)*y(2)+C*y(2)-A*y(1)*y(3);

y(5); (1/A1)*(A2*y(2)*y(4)-A3*y(1)*y(5)-A5*y(3)*y(3)-A4*y(4)-A6*y(2)*y(2))];

end

% Define the boundary conditions

function res = OdeBc(ya, yb)

res= [ya(1);ya(2)-1;ya(4)-1;yb(2);yb(4)];

end

% setting the initial guess for first solution

function v = OdeInit1(x,A)

v=[0.56;0;0;0;0];

end

[Torsten]

Use

sol = bvp4c (@(x,y)OdeBVP(x,y,A,B,C,A1,A2,A3,A4,A5,A6), @(ya,yb)OdeBc(ya, yb), solinit, options);
% Calculate the Nusselt number
Nu = -phi4*sqrt((aa+1/2)) * sol.y(5,:)
% Plot the Nusselt number
figure(3)
plot(sol.x, Nu, '-')
xlabel('\eta')
ylabel('Nu')
xlim([0 8])
ylim([0 1])
title('Nusselt Number Profile')

instead of

sol = bvp4c (@(x,y)OdeBVP(x,y,A,B,C,A1,A2,A3,A4,A5,A6), @(ya,yb)OdeBc(ya, yb), solinit, options);
eta = linspace (etaMin, etaMax1, stepsize1)
y= deval (sol,eta)
% Calculate the Nusselt number
Nu = -phi4*sqrt((aa+1/2)) * sol.y(5,:)
% Plot the Nusselt number
figure(3)
plot(sol.x, Nu, '-')
xlabel('\eta')
ylabel('Nu')
xlim([0 8])
ylim([0 1])
title('Nusselt Number Profile')

[Waseef]

Sir i want to plot this the nusselt number against the parameter (aa) which i mention in code but when i apply the loop it does not work how i do it

function waslipflow

format long g

%Define all parameters

% Boundary layer thickness & stepsize

etaMin = 0.01; etaMax1 = 10; etaMax2 = 15; stepsize1 = etaMax1;

stepsize2 = etaMax2;

% Input for the parameters

Phi=0.01;%input ('Input the value of phi = ');

R =0; %input ('Input the value of Radiation = ');

aa = 10;%input ('Input the value of aa = ');(Against this parameter)

M= 0.0;%input('magentic parameter M =');

pm=0.0;%input('porosity =');

Q=0.0; %input('Heat sourse parameter');

sigmaf=0.05;

sigmas=59600000;

n=3;%input ('Input the value of n = ');

ks= 400; kf= 0.613; Rhos= 8933; Rhof= 997.1;

Cps= 385; Cpf= 4179; tt=Rhos*Cps; kk=Rhof*Cpf;

Pr = 10; % input ('Input the Prandtl number = ');

Ec =0.0;%input ('Input Eckret for velociy exponent parameter = ');

spn=(3*Phi*((sigmas/sigmaf)-1))/(((sigmas/sigmaf)+2)-((sigmas/sigmaf)-1)*Phi);

sigmanf=sigmaf*(1+spn);

phi1=(1-Phi)^2.5;

phi2=1-Phi+Phi*(Rhos/Rhof);

A=phi1*phi2;

B=phi1*phi2*(2*aa/(aa+1));

C=(2/(aa+1))*((phi1*M*(sigmanf/sigmaf))+pm);

total=sigmas/sigmaf;

phi3=1+(3*(total-1)*Phi/((total+2)-(total-1)*Phi));

phi4=(ks+(n-1)*kf+(n-1)*(ks-kf)*Phi)/(ks+(n-1)*kf+(kf-ks)*Phi);

phi5=1-Phi+Phi*(Cps/Cpf);

phi5=1-Phi+Phi*(tt/kk);

A1=phi4+R;

A2=Pr*phi5*(4*aa/(aa+1));

A3=Pr*phi5;

A4=(2/(aa+1))*Pr*Q;

A5=Ec*Pr/phi1;

A6=(2/(aa+1))*Pr*Ec*M*(sigmanf/sigmaf);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%First solution %%%%%%%%%%%%%%%%%%%%%%

options = bvpset('stats','off','RelTol',1e-6);

solinit = bvpinit (linspace (etaMin, etaMax1, stepsize1),@(x)OdeInit1(x,A));

sol = bvp4c (@(x,y)OdeBVP(x,y,A,B,C,A1,A2,A3,A4,A5,A6), @(ya,yb)OdeBc(ya, yb), solinit, options);

% Calculate the Nusselt number

Nu = -phi4*sqrt((aa+1/2)) * sol.y(5,1);

% Plot the Nusselt number

figure(3)

plot(sol.x, Nu, '-')

xlabel('\eta')

ylabel('Nu')

xlim([0 8])

ylim([0 1])

title('Nusselt Number Profile')

end

% Define the ODE function

function f = OdeBVP(~,y,A,B,C,A1,A2,A3,A4,A5,A6)

f =[y(2); y(3);B*y(2)*y(2)+C*y(2)-A*y(1)*y(3);

y(5); (1/A1)*(A2*y(2)*y(4)-A3*y(1)*y(5)-A5*y(3)*y(3)-A4*y(4)-A6*y(2)*y(2))];

end

% Define the boundary conditions

function res = OdeBc(ya, yb)

res= [ya(1);ya(2)-1;ya(4)-1;yb(2);yb(4)];

end

% setting the initial guess for first solution

function v = OdeInit1(x,A)

v=[0.56;0;0;0;0];

end

[Torsten]

% Calculate the Nusselt number
Nu = -phi4*sqrt((aa+1/2)) * sol.y(5,:)

instead of

% Calculate the Nusselt number
Nu = -phi4*sqrt((aa+1/2)) * sol.y(5,1)

[Waseef]

This my code for generating skin fraction and nusselt number but i got both of an empty graphs please educate me why is this happening

slipflow()

function slipflow
format long g
%Define all parameters
% Boundary layer thickness & stepsize
global A Pr a b c C1 P1 C2 P2 K2 C3 P3 H1 H2 H3 Kh D1 beta
global D2 D3 D4 H4 H5 H6 H7 H8 B1 B2 B3 B S s1 s2 s3 lambda Re 
etaMin = 0;
etaMax1 = 15;
etaMax2 = 15; %15, 10
stepsize1 = etaMax1;
stepsize2 = etaMax2;
% Input for the parameters
Re=0.5;
A=0.6;  %velocity slip
B=0.2;  %thermal slip
beta=0.02;  %heat gen/abs
S=2.4;  %suction(2.3,2.4,2.5)
Pr=6.2;     %prandtl number
Lambda = -1.4:0.1:1;
for i = 1:numel(Lambda)  %stretching shrinking
    lambda = Lambda(i);
    a=0.01; %phil-1st nanoparticle concentration
    b=0.01; %(0.01,0.05)phi2-2nd nanoparticle concentration
    c=a+b;  %phi-hnf concentration of hybrid nanoparticle 
    %%%%%%%%%%% 1st nanoparticle properties (Al2O3)%%%%%%%%%%%%
    C1=765;
    P1=3970;
    K1=40;
    B1=0.85/((10)^5);
    s1=35*(10)^6; %MHD
    %%%%%%%%%%% 2nd nanoparticle properties (Cu)%%%%%%%%%%%%
    C2=385; %specific heat
    P2=8933; %density
    K2=400; %thermal conductivity
    B2=1.67/((10)^5);   %thermal expansion
    s2=(59.6)*(10)^6; %MHD
    %%%%%%%%%%% Base fluid properties %%%%%%%%%%%%
    C3=4179; %specific heat
    P3=997.1; %density
    K3=0.613; %thermal conductivity
    B3=21/((10)^5);   %thermal expansion
    s3=0.05; %MHD
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%multiplier%%%%%%%%%%%%%%%%%%%
    H1=P1*C1;       %pho*cp nanoparticle 1
    H2=P2*C2;       %pho*cp nanoparticle 2
    H3=P3*C3;       %pho*cp base fluid
    H4=a*H1+b*H2+(1-c)*H3; %pho*cp hybrid nanofluid
    H5=a*P1+b*P2+(1-c)*P3;  %pho hybrid nanofluid
    H6=1/((1-c)^2.5);   % mu hybrid nanofluid / mu base fluid
    H7=a*(P1*B1)+b*(P2*B2)+(1-c)*(P3*B3);   % thermal expansion of hybrid nanofluid
    %Kn=K3*(K1+2*K3-2*a*(K3-K1))/(K1+2*K3+a*(K3-K1));     %thermal conductivity of nanofluid
    Kh=(((a*K1+b*K2)/c)+2*K3+2*(a*K1+b*K2)-2*c*K3)/(((a*K1+b*K2)/c)+2*K3-(a*K1+b*K2)-2*c*K3); %khnf/kf
    H8=(((a*s1+b*s2)/c)+2*s3+2*(a*s1+b*s2)-2*c*s3)/(((a*s1+b*s2)/c)+2*s3-(a*s1+b*s2)-2*c*s3); % \sigma hnf/ \sigma f
    D1=(H5/P3)/H6;
    D3=(H7/(P3*B3))/(H5/P3);                % multiplier of boundary parameter
    D2= Pr*((H4/H3)/Kh);
    D4=H8/(H5/P3);  %multiplier MHD
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% First solution %%%%%%%%%%%%%%%%%%%
    options = bvpset('stats','off','RelTol',1e-10);
    solinit = bvpinit (linspace (etaMin, etaMax1, stepsize1),@(x)OdeInit1);
    sol=bvp4c(@OdeBVP, @OdeBC, solinit);
    To_Plot1(i) = H6*sol.y(3,1);
    To_Plot2(i) = Kh*sol.y(5,1);
    %sol = bvp4c (@(x,y)OdeBVP(x,y,Pr,D1,Kh,H4,H3,beta), @(ya,yb)OdeBc(ya, yb, A, S, B, lambda), solinit, options);
    %eta = linspace (etaMin, etaMax1, stepsize1);
    %y= deval (sol,eta);
    
    % saving the out put in text file for first solution
    descris =[sol.x; sol.y];
    %save 'sliphybrid_upper.txt' descris -ascii
    % Displaying the output for first solution
    % fprintf('\n First solution:\n');
    % fprintf('f"(0)=%7.9f\n',y(3));  % reduced skin friction 
    % fprintf('-theta(0)=%7.9f\n',-y(5)); %reduced local nusselt number
    % fprintf('Cfx=%7.9f\n',H6*(y(3)));    %  skin friction 
    % fprintf('Nux=%7.9f\n',-Kh*y(5));     % local nusselt number
    % fprintf('\n');
    %fprintf('%3.2f     %7.6f^\n',lambda,H6*y(3),-Kh*y(5))
end
figure(1)   %velocity profile
plot(Lambda,To_Plot1,'-')
xlabel('\lambda')
ylabel('f^\prime^\prime(0)')
xlim([-2.5 1])
figure(2)   %temperature profile
plot(Lambda,To_Plot2,'-')
xlabel('\lambda')
ylabel('\theta^\prime(0)')
xlim([-1.5 1])

% Define the ODE function
    function f = OdeBVP(~,y)
        f =[y(2);y(3);    D1*(2*(y(2)*y(2))-y(1)*y(3));y(5);(Pr/Kh)*((-H4/H3)*(y(1)*y(5)-y(2)*y(4))-beta*y(4))];
    end
% Define the boundary conditions
    function res = OdeBC(ya, yb)
        res= [ya(1)-S;ya(2)-lambda-A*ya(3);ya(4)-1-B*ya(5); yb(2); yb(4)];
    end
% setting the initial guess for first solution
    function v = OdeInit1(~)
        v=[S+0.56
            0
            0
            0
            0];
    end
end

[Torsten]

I corrected your code from above to make it work.

[Waseef]

Now its working thank you so much

[Waseef]

Hello, i want to geneate a plot of stream function in this code how i define it in this code

function slipflow

format long g

%Define all parameters

% Boundary layer thickness & stepsize

global A Pr a b c C1 P1 C2 P2 K2 C3 P3 H1 H2 H3 Kh D1 beta

global D2 D3 D4 H4 H5 H6 H7 H8 B1 B2 B3 B S s1 s2 s3 lambda Re

etaMin = 0;

etaMax1 = 15;

etaMax2 = 15; %15, 10

stepsize1 = etaMax1;

stepsize2 = etaMax2;

% Input for the parameters

Re=0.5;

A=0.6; %velocity slip

B=0.2; %thermal slip

beta=0.02; %heat gen/abs

S=2.4; %suction(2.3,2.4,2.5)

Pr=6.2; %prandtl number

Lambda = -1.4:0.1:1;

for i = 1:numel(Lambda) %stretching shrinking

lambda = Lambda(i);

a=0.01; %phil-1st nanoparticle concentration

b=0.01; %(0.01,0.05)phi2-2nd nanoparticle concentration

c=a+b; %phi-hnf concentration of hybrid nanoparticle

%%%%%%%%%%% 1st nanoparticle properties (Al2O3)%%%%%%%%%%%%

C1=765;

P1=3970;

K1=40;

B1=0.85/((10)^5);

s1=35*(10)^6; %MHD

%%%%%%%%%%% 2nd nanoparticle properties (Cu)%%%%%%%%%%%%

C2=385; %specific heat

P2=8933; %density

K2=400; %thermal conductivity

B2=1.67/((10)^5); %thermal expansion

s2=(59.6)*(10)^6; %MHD

%%%%%%%%%%% Base fluid properties %%%%%%%%%%%%

C3=4179; %specific heat

P3=997.1; %density

K3=0.613; %thermal conductivity

B3=21/((10)^5); %thermal expansion

s3=0.05; %MHD

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%multiplier%%%%%%%%%%%%%%%%%%%

H1=P1*C1; %pho*cp nanoparticle 1

H2=P2*C2; %pho*cp nanoparticle 2

H3=P3*C3; %pho*cp base fluid

H4=a*H1+b*H2+(1-c)*H3; %pho*cp hybrid nanofluid

H5=a*P1+b*P2+(1-c)*P3; %pho hybrid nanofluid

H6=1/((1-c)^2.5); % mu hybrid nanofluid / mu base fluid

H7=a*(P1*B1)+b*(P2*B2)+(1-c)*(P3*B3); % thermal expansion of hybrid nanofluid

%Kn=K3*(K1+2*K3-2*a*(K3-K1))/(K1+2*K3+a*(K3-K1)); %thermal conductivity of nanofluid

Kh=(((a*K1+b*K2)/c)+2*K3+2*(a*K1+b*K2)-2*c*K3)/(((a*K1+b*K2)/c)+2*K3-(a*K1+b*K2)-2*c*K3); %khnf/kf

H8=(((a*s1+b*s2)/c)+2*s3+2*(a*s1+b*s2)-2*c*s3)/(((a*s1+b*s2)/c)+2*s3-(a*s1+b*s2)-2*c*s3); % \sigma hnf/ \sigma f

D1=(H5/P3)/H6;

D3=(H7/(P3*B3))/(H5/P3); % multiplier of boundary parameter

D2= Pr*((H4/H3)/Kh);

D4=H8/(H5/P3); %multiplier MHD

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% First solution %%%%%%%%%%%%%%%%%%%

options = bvpset('stats','off','RelTol',1e-10);

solinit = bvpinit (linspace (etaMin, etaMax1, stepsize1),@(x)OdeInit1);

sol=bvp4c(@OdeBVP, @OdeBC, solinit);

To_Plot1(i) = H6*sol.y(3,1);

To_Plot2(i) = Kh*sol.y(5,1);

%sol = bvp4c (@(x,y)OdeBVP(x,y,Pr,D1,Kh,H4,H3,beta), @(ya,yb)OdeBc(ya, yb, A, S, B, lambda), solinit, options);

%eta = linspace (etaMin, etaMax1, stepsize1);

%y= deval (sol,eta);

% saving the out put in text file for first solution

descris =[sol.x; sol.y];

%save 'sliphybrid_upper.txt' descris -ascii

% Displaying the output for first solution

% fprintf('\n First solution:\n');

% fprintf('f"(0)=%7.9f\n',y(3)); % reduced skin friction

% fprintf('-theta(0)=%7.9f\n',-y(5)); %reduced local nusselt number

% fprintf('Cfx=%7.9f\n',H6*(y(3))); % skin friction

% fprintf('Nux=%7.9f\n',-Kh*y(5)); % local nusselt number

% fprintf('\n');

%fprintf('%3.2f %7.6f^\n',lambda,H6*y(3),-Kh*y(5))

end

figure(1) %velocity profile

plot(Lambda,To_Plot1,'-')

xlabel('\lambda')

ylabel('f^\prime^\prime(0)')

xlim([-2.5 1])

figure(2) %temperature profile

plot(Lambda,To_Plot2,'-')

xlabel('\lambda')

ylabel('\theta^\prime(0)')

xlim([-1.5 1])

% Define the ODE function

function f = OdeBVP(~,y)

f =[y(2);y(3); D1*(2*(y(2)*y(2))-y(1)*y(3));y(5);(Pr/Kh)*((-H4/H3)*(y(1)*y(5)-y(2)*y(4))-beta*y(4))];

end

% Define the boundary conditions

function res = OdeBC(ya, yb)

res= [ya(1)-S;ya(2)-lambda-A*ya(3);ya(4)-1-B*ya(5); yb(2); yb(4)];

end

% setting the initial guess for first solution

function v = OdeInit1(~)

v=[S+0.56

0

0

0

0];

end

end

[Torsten]

I don't know what you mean by

i want to geneate a plot of stream function in this code how i define it in this code

[Waseef]

sorry sir for inconveniance, how i generate stream lines for this problem like in the image thank you.

[Torsten]

In your code, η is between 0 and 15 instead of 0 and 1 and ξ is not defined.

Answer 1

[Torsten]

function slipflow
format long g
%Define all parameters
% Boundary layer thickness & stepsize
etaMin = 0;
etaMax1 = 15;
etaMax2 = 15; %15, 10
stepsize1 = etaMax1;
stepsize2 = etaMax2;
% Input for the parameters
A=0.6;  %velocity slip
B=0.2;  %thermal slip
beta=0.02;  %heat gen/abs
S=2.4;  %suction(2.3,2.4,2.5)
Pr=6.2;     %prandtl number
lambda=-1;  %stretching shrinking
a=0.01; %phil-1st nanoparticle concentration
b=0.01; %(0.01,0.05)phi2-2nd nanoparticle concentration
c=a+b;  %phi-hnf concentration of hybrid nanoparticle 
%%%%%%%%%%% 1st nanoparticle properties (Al2O3)%%%%%%%%%%%%
C1=765;
P1=3970;
K1=40;
B1=0.85/((10)^5);
s1=35*(10)^6; %MHD
%%%%%%%%%%% 2nd nanoparticle properties (Cu)%%%%%%%%%%%%
C2=385; %specific heat
P2=8933; %density
K2=400; %thermal conductivity
B2=1.67/((10)^5);   %thermal expansion
s2=(59.6)*(10)^6; %MHD
%%%%%%%%%%% Base fluid properties %%%%%%%%%%%%
C3=4179; %specific heat
P3=997.1; %density
K3=0.613; %thermal conductivity
B3=21/((10)^5);   %thermal expansion
s3=0.05; %MHD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%multiplier%%%%%%%%%%%%%%%%%%%
H1=P1*C1;       %pho*cp nanoparticle 1
H2=P2*C2;       %pho*cp nanoparticle 2
H3=P3*C3;       %pho*cp base fluid
H4=a*H1+b*H2+(1-c)*H3; %pho*cp hybrid nanofluid
H5=a*P1+b*P2+(1-c)*P3;  %pho hybrid nanofluid
H6=1/((1-c)^2.5);   % mu hybrid nanofluid / mu base fluid
H7=a*(P1*B1)+b*(P2*B2)+(1-c)*(P3*B3);   % thermal expansion of hybrid nanofluid
%Kn=K3*(K1+2*K3-2*a*(K3-K1))/(K1+2*K3+a*(K3-K1));     %thermal conductivity of nanofluid
Kh=(((a*K1+b*K2)/c)+2*K3+2*(a*K1+b*K2)-2*c*K3)/(((a*K1+b*K2)/c)+2*K3-(a*K1+b*K2)-2*c*K3); %khnf/kf
H8=(((a*s1+b*s2)/c)+2*s3+2*(a*s1+b*s2)-2*c*s3)/(((a*s1+b*s2)/c)+2*s3-(a*s1+b*s2)-2*c*s3); % \sigma hnf/ \sigma f
D1=(H5/P3)/H6;
D3=(H7/(P3*B3))/(H5/P3);                % multiplier of boundary parameter
D2= Pr*((H4/H3)/Kh);
D4=H8/(H5/P3);  %multiplier MHD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% First solution %%%%%%%%%%%%%%%%%%%
options = bvpset('stats','off','RelTol',1e-10);
solinit = bvpinit (linspace (etaMin, etaMax1, stepsize1),@(x)OdeInit1(x,A,S,lambda));
sol = bvp4c (@(x,y)OdeBVP(x,y,Pr,D1,Kh,H4,H3,beta), @(ya,yb)OdeBC(ya, yb, A, S, B, lambda), solinit, options);
eta = linspace (etaMin, etaMax1, stepsize1);
y= deval (sol,eta);
figure(1)   %velocity profile
plot(sol.x,sol.y(2,:),'-')
xlabel('\eta')
ylabel('f`(\eta)')
hold on
figure(2)   %temperature profile
plot(sol.x,sol.y(4,:),'-')
xlabel('\eta')
ylabel('\theta(\eta)')
hold on
% saving the out put in text file for first solution
descris =[sol.x; sol.y];
save 'sliphybrid_upper.txt' descris -ascii
% Displaying the output for first solution
fprintf('\n First solution:\n');
fprintf('f"(0)=%7.9f\n',y(3));  % reduced skin friction 
fprintf('-theta(0)=%7.9f\n',-y(5)); %reduced local nusselt number
fprintf('Cfx=%7.9f\n',H6*(y(3)));    %  skin friction 
fprintf('Nux=%7.9f\n',-Kh*y(5));     % local nusselt number
fprintf('\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% second solution %%%%%%%%%%%%%%%%%%%
options = bvpset('stats','off','RelTol',1e-10);
solinit = bvpinit (linspace (etaMin, etaMax2, stepsize2),@(x)OdeInit2(x,A,S,lambda));
sol = bvp4c (@(x,y)OdeBVP(x,y,Pr,D1,Kh,H4,H3,beta), @(ya,yb)OdeBC(ya, yb, A, S, B, lambda), solinit, options);
eta= linspace (etaMin, etaMax2, stepsize2);
y = deval (sol,eta);
figure(1)   %velocity profile
plot(sol.x,sol.y(2,:),'--')
xlabel('\eta')
ylabel('f`(\eta)')
hold on
figure(2)   %temperature profile
plot(sol.x,sol.y(4,:),'--')
xlabel('\eta')
ylabel('\theta(\eta)')
hold on
% saving the out put in text file for second solution
descris=[sol.x; sol.y];
save 'sliphybrid_lower.txt'descris -ascii
% Displaying the output for first solution
fprintf('\nSecond solution:\n');
fprintf('f"(0)=%7.9f\n',y(3));  % reduced skin friction 
fprintf('-theta(0)=%7.9f\n',-y(5)); %reduced local nusselt number
fprintf('Cfx=%7.9f\n',H6*(y(3)));    %  skin friction 
fprintf('Nux=%7.9f\n',-Kh*y(5));     % local nusselt number
fprintf('\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
% Define the ODE function
function f = OdeBVP(x,y,Pr,D1,Kh,H4,H3,beta)
f =[y(2);y(3);D1*(2*(y(2)*y(2))-y(1)*y(3));y(5);(Pr/Kh)*((-H4/H3)*(y(1)*y(5)-y(2)*y(4))-beta*y(4))];
end
% Define the boundary conditions
function res = OdeBc (ya, yb, A, S, B, lambda)
res= [ya(1)-S;ya(2)-lambda-A*ya(3);ya(4)-1-B*ya(5);yb(2);yb(4)];
end
% setting the initial guess for first solution
function v = OdeInit1(x,A,S,lambda)
v=[S+0.56;0;0;0;0];
end
% setting the initial guess for second solution
function v1 =OdeInit2(x, A, S,lambda)
v1 = [exp(-x);exp(-x);-exp(-x);-exp(-x);-exp(-x)];
end

[Torsten]

Rename the function

OdeBc

in

OdeBC

For MATLAB, small and capital letters define different things.

Same with the letter "s" you used in your original code. I hope you always meant your "S" - I replaced it where it was needed.

Answer 2

[Farooq Aamir]

This working now.

slipflow()

 First solution:
f"(0)=0.954395347
-theta(0)=3.658786212
Cfx=1.003836827
Nux=3.881552117

Second solution:
f"(0)=0.075669302
-theta(0)=3.617987101
Cfx=0.079589273
Nux=3.838268944

function slipflow
format long g
%Define all parameters
% Boundary layer thickness & stepsize
etaMin = 0;
etaMax1 = 15;
etaMax2 = 15; %15, 10
stepsize1 = etaMax1;
stepsize2 = etaMax2;
% Input for the parameters
A=0.6;  %velocity slip
B=0.2;  %thermal slip
beta=0.02;  %heat gen/abs
S=2.4;  %suction(2.3,2.4,2.5)
Pr=6.2;     %prandtl number
lambda=-1;  %stretching shrinking
a=0.01; %phil-1st nanoparticle concentration
b=0.01; %(0.01,0.05)phi2-2nd nanoparticle concentration
c=a+b;  %phi-hnf concentration of hybrid nanoparticle 
%%%%%%%%%%% 1st nanoparticle properties (Al2O3)%%%%%%%%%%%%
C1=765;
P1=3970;
K1=40;
B1=0.85/((10)^5);
s1=35*(10)^6; %MHD
%%%%%%%%%%% 2nd nanoparticle properties (Cu)%%%%%%%%%%%%
C2=385; %specific heat
P2=8933; %density
K2=400; %thermal conductivity
B2=1.67/((10)^5);   %thermal expansion
s2=(59.6)*(10)^6; %MHD
%%%%%%%%%%% Base fluid properties %%%%%%%%%%%%
C3=4179; %specific heat
P3=997.1; %density
K3=0.613; %thermal conductivity
B3=21/((10)^5);   %thermal expansion
s3=0.05; %MHD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%multiplier%%%%%%%%%%%%%%%%%%%
H1=P1*C1;       %pho*cp nanoparticle 1
H2=P2*C2;       %pho*cp nanoparticle 2
H3=P3*C3;       %pho*cp base fluid
H4=a*H1+b*H2+(1-c)*H3; %pho*cp hybrid nanofluid
H5=a*P1+b*P2+(1-c)*P3;  %pho hybrid nanofluid
H6=1/((1-c)^2.5);   % mu hybrid nanofluid / mu base fluid
H7=a*(P1*B1)+b*(P2*B2)+(1-c)*(P3*B3);   % thermal expansion of hybrid nanofluid
%Kn=K3*(K1+2*K3-2*a*(K3-K1))/(K1+2*K3+a*(K3-K1));     %thermal conductivity of nanofluid
Kh=(((a*K1+b*K2)/c)+2*K3+2*(a*K1+b*K2)-2*c*K3)/(((a*K1+b*K2)/c)+2*K3-(a*K1+b*K2)-2*c*K3); %khnf/kf
H8=(((a*s1+b*s2)/c)+2*s3+2*(a*s1+b*s2)-2*c*s3)/(((a*s1+b*s2)/c)+2*s3-(a*s1+b*s2)-2*c*s3); % \sigma hnf/ \sigma f
D1=(H5/P3)/H6;
D3=(H7/(P3*B3))/(H5/P3);                % multiplier of boundary parameter
D2= Pr*((H4/H3)/Kh);
D4=H8/(H5/P3);  %multiplier MHD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% First solution %%%%%%%%%%%%%%%%%%%
options = bvpset('stats','off','RelTol',1e-10);
solinit = bvpinit (linspace (etaMin, etaMax1, stepsize1),@(x)OdeInit1(x,A,S,lambda));
sol = bvp4c (@(x,y)OdeBVP(x,y,Pr,D1,Kh,H4,H3,beta), @(ya,yb)OdeBc(ya, yb, A, S, B, lambda), solinit, options);
eta = linspace (etaMin, etaMax1, stepsize1);
y= deval (sol,eta);
figure(1)   %velocity profile
plot(sol.x,sol.y(2,:),'-')
xlabel('\eta')
ylabel('f`(\eta)')
hold on
figure(2)   %temperature profile
plot(sol.x,sol.y(4,:),'-')
xlabel('\eta')
ylabel('\theta(\eta)')
hold on
% saving the out put in text file for first solution
descris =[sol.x; sol.y];
%save 'sliphybrid_upper.txt' descris -ascii
% Displaying the output for first solution
fprintf('\n First solution:\n');
fprintf('f"(0)=%7.9f\n',y(3));  % reduced skin friction 
fprintf('-theta(0)=%7.9f\n',-y(5)); %reduced local nusselt number
fprintf('Cfx=%7.9f\n',H6*(y(3)));    %  skin friction 
fprintf('Nux=%7.9f\n',-Kh*y(5));     % local nusselt number
fprintf('\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% second solution %%%%%%%%%%%%%%%%%%%
options = bvpset('stats','off','RelTol',1e-10);
solinit = bvpinit (linspace (etaMin, etaMax2, stepsize2),@(x)OdeInit2(x,A,S,lambda));
sol = bvp4c (@(x,y)OdeBVP(x,y,Pr,D1,Kh,H4,H3,beta), @(ya,yb)OdeBc(ya, yb, A, S, B, lambda), solinit, options);
eta= linspace (etaMin, etaMax2, stepsize2);
y = deval (sol,eta);
figure(1)   %velocity profile
plot(sol.x,sol.y(2,:),'--')
xlabel('\eta')
ylabel('f`(\eta)')
hold on
figure(2)   %temperature profile
plot(sol.x,sol.y(4,:),'--')
xlabel('\eta')
ylabel('\theta(\eta)')
hold on
% saving the out put in text file for second solution
descris=[sol.x; sol.y];
%save 'sliphybrid_lower.txt'descris -ascii
% Displaying the output for first solution
fprintf('\nSecond solution:\n');
fprintf('f"(0)=%7.9f\n',y(3));  % reduced skin friction 
fprintf('-theta(0)=%7.9f\n',-y(5)); %reduced local nusselt number
fprintf('Cfx=%7.9f\n',H6*(y(3)));    %  skin friction 
fprintf('Nux=%7.9f\n',-Kh*y(5));     % local nusselt number
fprintf('\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
% Define the ODE function
function f = OdeBVP(x,y,Pr,D1,Kh,H4,H3,beta)
f =[y(2);y(3);D1*(2*(y(2)*y(2))-y(1)*y(3));y(5);(Pr/Kh)*((-H4/H3)*(y(1)*y(5)-y(2)*y(4))-beta*y(4))];
end
% Define the boundary conditions
function res = OdeBc(ya, yb, A, S, B, lambda)
res= [ya(1)-S;ya(2)-lambda-A*ya(3);ya(4)-1-B*ya(5);yb(2);yb(4)];
end
% setting the initial guess for first solution
function v = OdeInit1(x,A,S,lambda)
v=[S+0.56;0;0;0;0];
end
% setting the initial guess for second solution
function v1 =OdeInit2(x, A, S,lambda)
v1 = [exp(-x);exp(-x);-exp(-x);-exp(-x);-exp(-x)];
end

[Nurain]

why i try to run this coding but still error at line 56

[Torsten]

Which MATLAB version do you use ?

What is line 56 ?

[khuram Rafique]

sol = bvp4c (@(x,y)OdeBVP(x,y,Pr,D1,Kh,H4,H3,beta), @(ya,yb)OdeBC(ya, yb, A, S, B, lambda), solinit, options);

[khuram Rafique]

still eror in line 56 which mentioned above

[Farooq Aamir]

Your error please?

[Farooq Aamir]

sol = bvp4c (@(x,y)OdeBVP(x,y,Pr,D1,Kh,H4,H3,beta), @(ya,yb)OdeBc(ya, yb, A, S, B, lambda), solinit, options); Use this it will work, kindly check your fuction name while calling it.

[khuram Rafique]

Thanks Farooq Aamir sir, its working now

[Ramanuja]

Thanks Farooq Aamir sir,

Thanks Farooq Aamir sir,

Thanks Farooq Aamir sir,

[Yasir]

Hello sir, What changes should i make in the code to plot the graphs of skin friction,Nusselt number and sherword number.

how can i get the different [oints to plot them

Answer 3

[Waseef]

sir how i define entropy in this code

this the entropy "NN=y(6)*y(6))+C1*y(7)+D*(y(5)*y(5)+y(2)*y(2))+E*(y(1)*y(8)+y(4)*y(8))"

and this is the code

Skinforbydirectional()

f"(0)=-1.145548435
g"(0)=-0.114554844
-theta(0)=1.029922860

function Skinforbydirectional
format long g
% Boundary layer thickness & stepsize
global A Pr aa pm Phi R M pm Q sigmaf sigmas sigmanf  n ks kf Rhos Rhof
global Cps Cpf tt kk Ec spn phi1 phi5 phi2 Lambda A B C st A2 A3 A4 A5 A6 total 
etaMin = 0;
etaMax = 10;
stepsize = etaMax;% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Define all parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Phi=0.04;          %input ('Input the value of phi = ');  
R =0.2;               %input ('Input the value of Radiation = '); 
M= 0.01;            %input('magentic parameter M =');
pm=0.2;             %input('porosity =');
Q=0.2;              %input('Heat sourse parameter');
%alpha=0.3;
aa=0.5;
%n=3;                %input ('Input the value of n = ');
Ec =0;        %input ('Input Eckret for velociy exponent parameter = ');
Pr = 6.8; % input ('Input the Prandtl number = ');
st=0.1;
%aa=10;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sigmaf=0.05; sigmas=1000000;
Ks= 400;   Kf= 0.613;
Rhos= 8933;  Rhof= 997.1;Cps= 385;  Cpf= 4179;  tt=Rhos*Cps;
kk=Rhof*Cpf;
spn=(3*Phi*((sigmas/sigmaf)-1))/(((sigmas/sigmaf)+2)-((sigmas/sigmaf)-1)*Phi);
sigmanf=sigmaf*(1+spn);
% Lambda = 1:1:10;
%for ii = 1:numel(Lambda)  %stretching shrinking
%     aa = Lambda(ii);
total=(sigmas/sigmaf);
phi3=1+(3*(total-1)*Phi/((total+2)-(total-1)*Phi));
%phi4=(ks+(n-1)*kf+(n-1)*(ks-kf)*Phi)/(ks+(n-1)*kf+(kf-ks)*Phi);
phi4=((1-Phi)+2*Phi*Ks/(Ks-Kf)*log((Ks+Kf)/2*Kf))/((1-Phi)+2*Phi*Kf/(Ks-Kf)*log((Ks+Kf)/2*Kf));
%phi5=1-Phi+Phi*(Cps/Cpf);
phi1=(1-Phi)^2.5;
phi2=1-Phi+Phi*(Rhos/Rhof);
phi5=1-Phi+Phi*(tt/kk);
A = ((phi1 * M * (sigmanf / sigmaf)) + pm) * (2 / (aa + 1));
B = phi1 * phi2 * (2 * aa / (aa + 1));
C = phi1 * phi2;
A1 = phi4 + R;
B1 = Ec * Pr / phi1;
C1 = Pr * Q * (2 / (aa + 1));
E = Pr * phi5;
D = Pr * Ec * M * (sigmanf / sigmaf) * (2 / (aa + 1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% First solution %%%%%%%%%%%%%%%%%%%
options = bvpset('stats','off','RelTol',1e-10);
solinit = bvpinit (linspace (etaMin, etaMax, stepsize),@(x)OdeInit1);
sol=bvp4c(@OdeBVP, @OdeBC, solinit);
%sol = bvp5c (@(x,y)OdeBVP(x,y,Pr,D1,Kh,H4,H3,beta), @(ya,yb)OdeBc(ya, yb, A, S, B, lambda), solinit, options);
eta = linspace (etaMin, etaMax, stepsize);
y= deval (sol,eta);
% To_Plot1(ii) = (1/phi1)*(sqrt(2*(aa+1)))*sol.y(5,1);
% To_Plot2(ii) = -(phi4+R)*(sqrt((aa+1)/2))*sol.y(8,1);
% fprintf('y(3) at etaMax = %f\n', y(3, end));
%     fprintf('y(8) at etaMax = %f\n', y(8, end));
%end
figure(1)   %velocity profile
plot(sol.x,sol.y(2,:),'-')
xlabel('\eta')
ylabel('f`(\eta)')
hold on
figure(2)   %velocity profile
plot(sol.x,sol.y(5,:),'-')
xlabel('\eta')
ylabel('q`(\eta)')
hold on
% figure(3)   %velocity profile
% plot(Lambda,To_Plot1,'LineWidth',2)
% xlabel('a')
% ylabel('f^\prime^\prime(0)')
% grid on
% hold on
% figure(4)   %temperature profile
% hold on
% grid on
% plot(Lambda,To_Plot2,'LineWidth',2)
% xlabel('a')
% ylabel('\theta^\prime(0)')
%  %xlim([0 2])
% figure(1)   %velocity profile
% plot(Lambda,To_Plot1,'-')
% xlabel('\lambda')
% ylabel('f^\prime^\prime(0)')
% xlim([1,2])
% figure(2)   %temperature profile
% plot(Lambda,To_Plot2,'-')
% xlabel('\lambda')
% ylabel('\theta^\prime(0)')
% xlim([1 2])
% Define the ODE function
fprintf('f"(0)=%7.9f\n',y(3));  % reduced skin friction 
fprintf('g"(0)=%7.9f\n',y(6));  % reduced skin friction 
fprintf('-theta(0)=%7.9f\n',-y(8)); %reduced local nusselt number
    function f = OdeBVP(~,y)
        f =[ y(2); y(3);A*y(2)+B*(y(2)*y(2)+y(5)*y(2))-C*(y(1)*y(3)+y(4)*y(3));
            y(5);y(6); A*y(5)+B*(y(5)*y(2)+y(5)*y(5))-C*(y(1)*y(6)+y(4)*y(6));
            y(8);(-1/A1)*(B1*(y(3)*y(3)+y(6)*y(6))+C1*y(7)+D*(y(5)*y(5)+y(2)*y(2))+E*(y(1)*y(8)+y(4)*y(8)))];
    end
% Define the boundary conditions
    function res = OdeBC(ya, yb)
        res= [ya(1);
            ya(2)-1;
            ya(4);
            ya(5)-st;
            ya(7)-1; 
            yb(2); 
            yb(5);
            yb(7)];
    end
% setting the initial guess for first solution
    function v = OdeInit1(~)
        v=[0.9
            0.1
            0.1
            -0.1
            0.1
            0
            -0.1
            01];
    end
end