how to solve this bvp4c?

Hello, can anyone help me solve this bvp4c program?

as to solve this problem, I already try numerous initial numbers.

it's either

Unable to solve the collocation equations -- a singular Jacobian encountered.

or

Warning: Unable to meet the tolerance without using more than 1250 mesh points.
 The last mesh of 892 points and the solution are available in the output argument.
 The maximum residual is 85.3834, while requested accuracy is 1e-07. 

here is the bvp4c:

format long g
global alpha beta phi phid epsilon delta rhoS rhoF Y betaT Ec cpS cpF kf knf Pr Tv A B w Q ks

a=0; b=15;
phi=0.5;phid=0.2;delta=0.5;A=0.5;B=0.5;w=1.5;
betaT=0.2;alpha=0.5;beta=0.5;Ec=0.2;Y=0.01;Pr=6.2;Tv=0.488;epsilon=0.1; 
rhoF=997.1;cpF=4179;kf=0.613;   
rhoS=8933;cpS=385;knf=0.81629325;  
Q=((ks+(2*kf))-(2*phi*(kf-ks)))/((ks+(2*kf))+(phi*(kf-ks))); ks=400;

solinit = bvpinit(linspace(a,b,100),@fluidparticle_init);
options = bvpset('stat','on','RelTol',1e-7); 
sol = bvp4c(@fluidparticle_ode,@fluidparticle_bc,solinit,options);

sol.y(3,1)      
sol.y(7,1)      
plot(sol.x,sol.y(2,:),':r')    

function dydx=fluidparticle_ode(x,y,alpha,beta,phi,phid,Y,betaT,Ec,kf,knf,Pr,Tv,A,B,rhoS,rhoF,cpS,cpF,Q)
global alpha beta phi phid Y betaT Ec kf knf Pr Tv A B rhoS rhoF cpS cpF Q
dydx=[y(2);
      y(3);
      (((1-phi)^2.5)*(1-phi+(phi*(rhoS/rhoF)))*((y(2)^2)-(y(1)*y(3))))-((((1-phi)^2.5)/(1-phid))*alpha*beta*(y(5)-y(2)));
      y(5);
      ((y(5)^2)-(beta*(y(5)-y(2))))/y(4);
      y(7);
        (1-phi+(phi*((rhoS*cpS)/(rhoF*cpF))))*(Pr/Q)*((2*y(2)*y(6))-(y(1)*y(7)))-((Pr/Q)*((alpha*betaT*(y(8)-y(6)))-((alpha*Ec/Tv)*((y(5)-y(2))^2))))
      ((2*y(5)*y(8))+(Y*betaT*(y(8)-y(6))))/y(4)];
function res=fluidparticle_bc(ya,yb,epsilon,delta,w)
global epsilon delta w
res=[ya(2)-epsilon-(delta*ya(3))
    ya(1)
    yb(2)
    yb(5)
    yb(4)-yb(1)
    ya(6)-1-(w*ya(7))
    yb(6)
    yb(8)];
function v=fluidparticle_init(x)
v=[0.5
   -exp(-x) 
  exp(-x)
   0.3
  -0.5
   -exp(-x)
   exp(-x)
   -exp(-x)];

what should I adjust in this program?

thank you in advance for your help and kindness.

[darova]

Please attach the original equations

[NURUL AISYAH JOHAN]

〖(1-ϕ_d)/(1-ϕ)^2.5〗f^''' - (1-ϕ_d) [1-ϕ+ϕ(ρ_s/ρ_f )](f^'2-〖ff〗^'' )+αβ(F^'- f^' )= 0.....(1)

F^'2-FF^''-β(f^'-F^' )=0.....(2)

(1/Pr)(k_nf/k_f ) θ^''-(1-ϕ+ϕ((ρC_p )_(s )/(ρC_p )_f ))(2f^' θ- fθ^' )+ αβ_T (θ_p-θ)+αEc/τ_v (F^'-f^' )^2 = 0.....(3)

2F^' θ_p-Fθ_p^'+γβ_T (θ_p-θ)=0.....(4)

Boundary condition:

f^' (0) = ε+δf^'' (0), f(0) = 0, θ(0) = 1+ωθ'(0)

f^' (η)=0, F^' (η)=0, F(η)=f(η), θ(η)=0, θ_p (η)=0 as η → ∞

Answer 1

[darova]

I re-wrote the main equations (compare yours and mine) and tried b = 1

Check the attachment

[NURUL AISYAH JOHAN]

Thank you so much for your help. The program can run but may I ask. Why the limit is only at b=1? the limit seem too small to analyse the graph. Can anything to be done about that?

[darova]

• Why the limit is only at b=1?

I tried b=2 and this is what i got

[NURUL AISYAH JOHAN]

thank you so much for your help. you have help me big time.

but, i have new question about this matter. Why I get error when apply negative value for epsilon?

when epsilon=-0.5, the result:

Error using bvp4c (line 251)

Unable to solve the collocation equations -- a singular Jacobian encountered.

Error in f (line 26)

sol = bvp4c(@fluidparticle_ode,@fluidparticle_bc,solinit,options);

[darova]

It's hard to say without knowing of the equations. What is the meaning of epsilon? Maybe it can't be negative?

[Didarul Ahasan Redwan]

why I wasnt able to set upper limit b=6/7?

thanks in advance

[darova]

Why wasn't you? I don't know

I set up a limit b=6/7 and got this

[Didarul Ahasan Redwan]

https://uk.mathworks.com/matlabcentral/answers/516147-boundary-value-problem-with-bvp4c

can you please have a look at this?

I would be very thankful