how can i plot graph of skin friction using the bellow code.

what changes are required in this code

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=1; %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.9fn’,y(3)); % reduced skin friction
fprintf(‘-theta(0)=%7.9fn’,-y(5)); %reduced local nusselt number
fprintf(‘Cfx=%7.9fn’,H6*(y(3))); % skin friction
fprintf(‘Nux=%7.9fn’,-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.9fn’,y(3)); % reduced skin friction
fprintf(‘-theta(0)=%7.9fn’,-y(5)); %reduced local nusselt number
fprintf(‘Cfx=%7.9fn’,H6*(y(3))); % skin friction
fprintf(‘Nux=%7.9fn’,-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)];
endwhat changes are required in this code

Cart

Cart