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

[Walter Roberson]

No figure was attached

[MINATI]

@ Walter

Sorry for the mistake (I have tried but it was not attached)

Now again I have tried so many times to attach the paper but unable.

I think you dont mind downloading the paper from sci-hub.tw removing barriers

with

doi.org/10.1155/2013/724547

[Jan]

@MINATI: It is not hard to attach a file. Please mention, which problem you have with it. I've attached it here now, but as soon as other comments are posted, it will run out of view. So please try again to attach it to the original question, where readers expect it: Open the question for editing, press the paperclip in the toolbar, select the file, press "Open" in the dialog, an click on "Submit" finally.

What is the problem with the posted code?

[MINATI]

Actually the attach symbol (clip symbol) is deactivated perhaps

[MINATI]

@ Jan I need to draw fig. 1

Actually again I repeat the same procedure but unable to attach.

Sorry again.

[Jan]

@MINATI: "I'm unable to attach" does not explain, what you are doing and which error message you get or which other problem occurs.

As far as I remember the clip symbol is disabled, if you have posted a certain number of attachments at a specific day already. You can simply find out if it is "deactivated perhaps" by clicking on it. What do you see then? (It is the one left beside the question mark.)

I've added the paper to the question now.

[MINATI]

Thanks Jan for the help

But see this time here it is attached

Any way

Please somebody follow the code and Fig. 1

[Walter Roberson]

Dig out the relevant equations and post them. Indicate the initial conditions. Indicate whether you are looking for numeric solution or theoretical solution. Figure out whether it is a "stiff" system, or if you need a Mass Matrix. I see an f''' in there: does that hint that we could provide the Jacobian or Hessian as functions?

Do you have the Symbolic Toolbox?

[MINATI]

@Walter

Relevant equations are in (7) and boundary conditions are in (8).

Yes we need numeric solution.

But can it be possible to derive theoretical solution from the code.

whether it is a "stiff" system. I dont know what is "stiff".

does that hint that we could provide the Jacobian or Hessian as functions?

Ans: I think its NO

[David Wilson]

If I understand correctly, the BVP you are trying to solve has BCs at infinity. You have chosen 6 (& the paper uses 8), so you might like to validate that approximation.

What exactly is the problem ?

My solution is for the Pr=6.2, \gamma=0.1. This seems to follow your Fig 1.

[Jan]

@MINATI: See https://en.wikipedia.org/wiki/Stiff_equation

Answer 1

[David Wilson]

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.

[MINATI]

Many many thanks David

It worked

Where to put the loop for Gamma G=[0 0.1 0.2 0.3] which will vary the fig