I'm trying too solve a set of coupled ODEs using BVP5C but I keep getting the same (wrong) results for the different cases. What am I doing wrong?

I'm trying to solve a set of coupled ODEs using BVP4C but the results I get from Matlab isn't making sense. I keep getting the same value.

n  = 0.4;
Bo = [0.01322917839
0.01082513952
0.004466835920
0.001258925411
0.0003548133892];
dy    = zeros(5,1);
dy(1) = y(2);
dy(2) = ((2*n+1)*y(3)*y(2)*Bo^(2/(n+1))/(n+1))/0.2e1;
dy(3) = y(4);
dy(4) = y(5);
dy(5) = ((-y(4)^2*n+2*y(3)*y(5)*n+2*y(1)*n-y(4)^2+y(3)*y(5)+2*y(1))/y(5)^(n-1)/n/(n+1))/0.2e1;

Where,

y(1) = theta

y(2) = theta'

y(3) = f

y(4) = f'; and

y(5) = f''

Subject to the following Boundary Conditions,

3 initial conditions are given: eta=0, f(0)= f'(0)=0,theta(0)=1

The boundary conditions that needs to be satisfied are: f'(eta=inf)= 0 and theta(eta=inf)=0.

I have tried the following code but I get f''= 0 and theta'= -0.033333333333333 for all cases.

    function sol= proj
clc;clf;clear;
global Bo n
% Bo=1;
% n=1;
pp=[0.4];
qq=[0.01322917839
0.01082513952
0.004466835920
0.001258925411
0.0003548133892];
%figure(1)
%plot(2,1);hold on
options=bvpset('stats','on','RelTol',1e-9);
m=linspace(0,30);
solinit= bvpinit(m,[1,0,0,0,0]);
for i=1:numel(pp)
    n=pp(i);
    for i=1:numel(qq);
        Bo=qq(i)
sol= bvp5c(@projfun,@projbc,solinit,options);
format long
 y1=deval(sol,0)
solinit= sol;
plot(sol.x,sol.y(1,:));hold on
end
end
end
function f= projfun(x,y)
global Bo n
f    = zeros(5,1);
f    = zeros(5,1);
f(1) = y(2);
f(2) = ((2*n+1)*y(3)*y(2)*Bo^(2/(n+1))/(n+1))/0.2e1;
f(3) = y(4);
f(4) = y(5);
f(5) = ((-y(4)^2*n+2*y(3)*y(5)*n+2*y(1)*n-y(4)^2+y(3)*y(5)+2*y(1))/y(5)^(n-1)/n/(n+1))/0.2e1;
end
function res= projbc(ya,yb)
res= [ya(1)-1; ya(3); ya(4); yb(1); yb(4)];
end

The code works perfectly for the following conditions when n =1 but gives the same wrong values for n other than 1.

n  = 1.0;
Bo = [0.75
0.5
0.1
0.01
0.001
];
dy    = zeros(5,1);
dy(1) = y(2);
dy(2) = ((2*n+1)*y(3)*y(2)*Bo^(2/(n+1))/(n+1))/0.2e1;
dy(3) = y(4);
dy(4) = y(5);
dy(5) = ((-y(4)^2*n+2*y(3)*y(5)*n+2*y(1)*n-y(4)^2+y(3)*y(5)+2*y(1))/y(5)^(n-1)/n/(n+1))/0.2e1;

What I'm I doing wrong?

[darova]

[Mayokun Ojediran]

@darova the value for y(5) isn't meant to be zero, but the code keeps giving zero.

[Torsten]

You prescribe y(5)=0 in the line

solinit= bvpinit(m,[1,0,0,0,0]);

[Mayokun Ojediran]

That's just an initial guess

[Torsten]

And what do you think BVP does with this guess ?

It starts with y(1)=1, y(2)=0, y(3)=0, y(4)=0 and y(5)=0, calls "projfun" and - oops - divides by zero.

[Mayokun Ojediran]

I have tried other values for the initial guess, still having the same problem.

[Torsten]

Use

 for j=1:numel(qq);
    Bo=qq(j)

instead of

 for i=1:numel(qq);
    Bo=qq(i)

[Mayokun Ojediran]

Done that, no change.

[darova]

Is this correct form for initial condition?

[Mayokun Ojediran]

Solinit is a function name in the BVP4C command