MATLAB: How to draw Fig. 1 from the attached pdf with this code

how to draw fig. 1 from the attached pdf with this code

function main
Pr=1; G=0.1;
% phi=input('phi=');   %%0,.05, .1, .15, .2
phi=0.0;
rhof=997.1;Cpf=4179;kf=0.613;           %for    WATER
rhos=6320;Cps=531.8;ks=76.5;               %for    CuO
a1=((1-phi)^2.5)*(1-phi+phi*(rhos/rhof));
a2=(1-phi+phi*((rhos*Cps)/(rhof*Cpf)));
A=(ks+2*kf+phi*(kf-ks))/(ks+2*kf-2*phi*(kf-ks));  %%%%Knf
xa=0;xb=6;
solinit=bvpinit(linspace(xa,xb,101),[0 1 0 1 0]);
sol=bvp4c(@ode,@bc,solinit);
xint=linspace(xa,xb,101);
sxint=deval(sol,xint);
figure(1)
plot(xint,(1-phi)^-2.5*sxint(3,:),'-','Linewidth',1.5);   %for f''(0)/(1-phi)^2.5  vs  phi
xlabel('\eta'); 
ylabel('f''(0)/(1-phi)^2.5');
hold on
function res=bc(ya,yb)
res=[ya(1); ya(2)-1-G*ya(3); ya(4)-1; yb(2); yb(4)];  
end
 function dydx=ode(x,y)    
dydx=[y(2);  y(3);  a1*(y(2)^2-y(3)*y(1));  y(5);  -A*Pr*a2*y(1)*y(5)];              
end
end
http://doi.org/10.1155/2013/724547
[EDITED, Jan, Attachment added].

Best Answer

I didn't bother draw the other 3 lines, but you just ned to make the necessary changes to gamma for that.
If you run something like what you had originally, you only want the fist point of f''(). 
Pr=6.2; G=0.1;
% phi=input('phi=');   %%0,.05, .1, .15, .2
phi=0.0;
rhof=997.1;Cpf=4179;kf=0.613;           %for    WATER
rhos=6320;Cps=531.8;ks=76.5;               %for    CuO
a1=((1-phi)^2.5)*(1-phi+phi*(rhos/rhof));
a2=(1-phi+phi*((rhos*Cps)/(rhof*Cpf)));
A=(ks+2*kf+phi*(kf-ks))/(ks+2*kf-2*phi*(kf-ks));  %%%%Knf

BCres= @(ya,yb) ... 
       [ya(1); ya(2)-1-G*ya(3); ya(4)-1; yb(2); yb(4)];  
fODE = @(x,y) ... 
       [y(2);  y(3);  a1*(y(2)^2-y(3)*y(1));  y(5);  -A*Pr*a2*y(1)*y(5)];  
xa=0;xb=8;
solinit=bvpinit(linspace(xa,xb,101),[0 1 0 1 0]);
sol=bvp4c(fODE,BCres,solinit);
xint=linspace(xa,xb,101);
sxint=deval(sol,xint);
figure(1)
plot(xint,(1-phi)^-2.5*sxint(3,:),'-','Linewidth',1.5);   %for f''(0)/(1-phi)^2.5  vs  phi
xlabel('\eta'); 
ylabel('f''(0)/(1-phi)^2.5');
Now you have to re-run the above, but change phi over the range given in the Fig. 
xa=0;xb=8;
phiv = [0:0.04:0.2]'; 
p = []; % collect points here 
for i=1:length(phiv)
    phi = phiv(i); 
    a1=((1-phi)^2.5)*(1-phi+phi*(rhos/rhof));
    a2=(1-phi+phi*((rhos*Cps)/(rhof*Cpf)));
    A=(ks+2*kf+phi*(kf-ks))/(ks+2*kf-2*phi*(kf-ks));  %%%%Knf
    
    fODE = @(x,y) ... 
       [y(2);  y(3);  a1*(y(2)^2-y(3)*y(1));  y(5);  -A*Pr*a2*y(1)*y(5)];  
    solinit=bvpinit(linspace(xa,xb,101),[0 1 0 1 0]);
    sol=bvp4c(fODE,BCres,solinit);
    
    p(i,1) = (1-phi)^-2.5*sxint(3,1)
end 
plot(phiv, p,'o-')
xlabel('\phi'); ylabel('f''''(0) & stuff')
Resultant plot is as above.
Related Question