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

what changes are required in this code

slipflow()

 First solution:
f"(0)=0.699647204
-theta(0)=3.665712766
Cfx=0.735891715
Nux=3.888900396

Warning: Attempt to write an unsupported data type to an ASCII file.
Variable 'options' not written to file.

Warning: Attempt to write an unsupported data type to an ASCII file.
Variable 'sol' not written to file.

Warning: Attempt to write an unsupported data type to an ASCII file.
Variable 'solinit' not written to file.

Second solution:
f"(0)=0.053709708
-theta(0)=3.618502237
Cfx=0.056492084
Nux=3.838815445

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.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

Answer 1

[Image Analyst]

I always use the function version of save. Try

save('sliphybrid_upper.txt', 'descris',  '-ascii');

[Yasir]

Then how to plot from the file