Shooting Method Boundary Conditions not Implementing

I am trying to implement RK4 using the shooting method for the following boundary layer equations,

n  = 0.4;
Bo = 0.01322917839;
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 as eta=30.

I have tried the following code but the graph I get is not as it's supposed to be.

function Shooting_Method_code1
clc
clear all
x = [0.5 0.5];
options = optimset('Display', 'iter');
x1 = fsolve(@solver,x);
end

function F = solver(x)
options = odeset ('RelTol', 1e-8, 'AbsTol', [1e-8 1e-8 1e-8 1e-8 1e-8]);
[t,u] = ode45(@equation, [0,1], [1 x(1) 0 0 x(2)], options);
s = length(t);
F = [u(s,1)-0, u(s,4)-0];

figure(1)
plot(t,u(:,1), t, u(:,4))
end

function dy = equation(t,y)
dy    = zeros(5,1);
n  = 0.4;
Bo = 0.01322917839;
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;
end

Please help.

[Torsten]

You provide f'(0)=0 while in your mathematical description you claim f'(0)=1.

[Mayokun Ojediran]

Torsten, that was an errror and has been corrected.

[Torsten]

The boundary conditions seem to be satisfied for the solution you get.

Why do you write that the graph is not as it is supposed to be ?

[Mayokun Ojediran]

The graphs should look something like this. I think the problem is the span from of the graph, but increasing tspan causes the code to fail.

[Torsten]

Did you try bvp4c ?

[Mayokun Ojediran]

I am not too familiar with bvp4c, could you please help out with it?

[Torsten]

function main
xmesh = linspace(0, 30, 1000);
solinit = bvpinit(xmesh, [0; 0; 0; 0; 0]);
sol = bvp4c(@equation, @bcfcn, solinit);
plot(sol.x,sol.y(4,:))
end

function dy = equation(t,y)
dy    = zeros(5,1);
n  = 0.4;
Bo = 0.01322917839;
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;
end

function res = bcfcn(ya,yb)
res(1) = ya(1)-1.0;
res(2) = ya(3);
res(3) = ya(4);
res(4) = yb(1);
res(5) = yb(4);
end

[Mayokun Ojediran]

Thank you, this works.