Numerically solving a 6th order non-linear differential equation in Matlab

Question

I've posted yesterday a question about solving a non linear equation : it was not clear so I am reformulating my question.

I am trying to solve a high-order non linear differential equation presented right below with $F \in [-\eta_{0}, \eta_{0}]$

$\left(5F - 4\eta\frac{dF}{d\eta}\right) = 9 \frac{d}{d\eta}\left(F^{3} \frac{d^{5}F}{d\eta^{5}}\right)$

As boundary conditions, I have the following ones :

$\left\{ \begin{array}{ll} F' = F'''=0 & \mbox{in } \space \eta = 0 \\ F = F' = 0 & \mbox{in} \space \eta =\pm \space\eta_{0} \end{array} \right.$

I am trying to apply a shooting method algorithm in order to solve this equation on $[0, \eta_{0}]$, I have converted it in a system of 1-st order differential equations as following :

$f_{i}' = \left[\begin{array}{@{}c@{}} y_{2}=\frac{dF}{d\eta} \\ y_{3} = \frac{d^{2}F}{d\eta^{2}}\\ y_{4} = \frac{d^{3}F}{d\eta^{3}} \\ y_{5} = \frac{d^{4}F}{d\eta^{4}} \\ y_{6} = \frac{d^{5}F}{d\eta^{5}}\\ y_{7} = \frac{d^{6}F}{d\eta^{6}} \end{array} \right]$ with : $y_{7} = \left[\left(\frac{5}{9y_{1}^{2}}\right) - \left(\frac{4xy_{2}}{9y_{1}^{3}}\right)-\left(\frac{3y_{2}y_{6}}{y_{1}}\right)\right]$

According to the boundary conditions we thus have : $y_{2}(0) = y_{4}(0) = 0$ and we have to guess the remaining ones : $y_{1}(0), y_{3}(0),y_{5}(0) \space \mbox{and} \space y_{6}(0)$

I have attached to my post, the algorithms that I wrote which are inspired by biblographical researchs : the main code is called "shoot4nl.m" and I have tried to guess some initial values in order to run the algorithm. Also, I have choses a value of $\eta_{0}$ very "high" in order to match the conditions...

Here are the algorithms that I am using : These are MATLAB scripts so I am sharing it below. The main function is the one called "shoot4nl".

function shoot4nl

% Shooting method for nonlinear 4th-order boundary
% value problem 

global XSTART XSTOP H   % Make these params. global.
XSTART = 0; XSTOP = 50; % Range of integration.
H = 0.25;                % Step size.
freq = 1;               % Frequency of printout.
u = [1 -1 1 0];             % Trial values of u(1)
                        % and u(2), u(3) and u(4).

x = XSTART;
u = newtonRaphson2(@residual,u);
[xSol,ySol] = runKut5(@dEqs,x,inCond(u),XSTOP,H);
printSol(xSol,ySol,freq)

function F = dEqs(x,y)    % Differential equations.
F = zeros(1,6);
F(1) = y(2); F(2) = y(3); F(3) = y(4); F(4) = y(5); F(5) = y(6);
F(6) = ((5/(9*y(1)^2))-(4*x*y(2)/(9*y(1)^3))-(3*y(2)*y(6)/y(1)));



function y = inCond(u)    % Initial conditions; u(1),u(2),u(3) and u(4) are unknowns
y = [u(1) 0 u(2) 0 u(3) u(4)];      


function r = residual(u)  % Bounday residuals.
global XSTART XSTOP H
r = zeros(length(u),1);
x = XSTART;
[xSol,ySol] = runKut5(@dEqs,x,inCond(u),XSTOP,H);
lastRow = size(ySol,1);
r(1) = ySol(lastRow,3);
r(2) = ySol(lastRow,4) - 1;

In order to run it we have to define other functions.

runKut5 :

function [xSol,ySol] = runKut5(dEqs,x,y,xStop,h,eTol)

% 5th-order Runge-Kutta integration.
% USAGE: [xSol,ySol] = runKut5(dEqs,x,y,xStop,h,eTol)
% INPUT:
% dEqs  = handle of function that specifyies the 1-st order differential
% equations : F(x,y) = [dy1/dx dy2/dx dy3/dx ...]

%
% x,y
% xStop = terminal value of x.
% h     = trial value of increment of x.
% eTol  = per-step error tolerance (default = 1.0e-6).
% OUTPUT:
% xSol = x-values at which solution is computed.
% ySol = values of y corresponding to the x-values.

if size(y,1) > 1 ; y = y'; end % y must be row vector

if nargin < 6; eTol = 1.0e-6; end

n = length(y);
A = [0 1/5 3/10 3/5 1 7/8]; 
B=[    0          0       0          0           0
       1/5        0       0          0           0
       3/40      9/40     0          0           0
       3/10     -9/10    6/5         0           0
     -11/54      5/2   -70/27      35/27         0
    1631/55296 175/512 575/13824 44275/110592 253/4096];
C = [37/378 0 250/621 125/594 0 512/1771];
D = [2825/27648 0 18575/48384 13525/55296 277/14336 1/4];
% Initialize solution
xSol = zeros(2,1); ySol = zeros(2,n);
xSol(1) = x; ySol(1,:) = y;
stopper = 0; k = 1;
for p = 2:5000
    % Compute K’s from Eq. (7.18)
    K = zeros(6,n);
    K(1,:) = h*feval(dEqs,x,y);
    for i = 2:6
        BK = zeros(1,n);
        for j = 1:i-1
            BK = BK + B(i,j)*K(j,:);
        end
        K(i,:) = h*feval(dEqs, x + A(i)*h, y + BK);
    end
    % Compute change in y and per-step error from
    % Eqs.(7.19) & (7.20)
    dy = zeros(1,n); E = zeros(1,n);
    for i = 1:6
        dy = dy + C(i)*K(i,:);
        E = E + (C(i) - D(i))*K(i,:);
    end
    e = sqrt(sum(E.*E)/n);
    % If error within tolerance, accept results and
    % check for termination
    if e <= eTol
        y = y + dy; x = x + h;
        k = k + 1;
        xSol(k) = x; ySol(k,:) = y;
        if stopper == 1;
            break 
        end
    end
    % Size of next integration step from Eq. (7.24)
    if e~=0; hNext = 0.9*h*(eTol/e)^0.2;
    else; hNext=h;
    end
    % Check if next step is the last one (works
    % with positive and negative h)
    if (h > 0) == (x + hNext >= xStop )
        hNext = xStop - x; stopper = 1;
    end
h = hNext; 
end

newtonRaphson2 :

function root = newtonRaphson2(func,x,tol)

% Newton-Raphson method of finding a root of simultaneous
% equations fi(x1,x2,...,xn) = 0, i = 1,2,...,n.
% USAGE: root = newtonRaphson2(func,x,tol)
% INPUT:
% func = handle of function that returns[f1,f2,...,fn].
% x    = starting solution vector [x1,x2,...,xn].
% tol  = error tolerance (default is 1.0e4*eps).
% OUTPUT:
% root = solution vector.

if nargin == 2
    tol = 1.0e4*eps; 
end

if size(x,1) == 1; x = x'; end   % x must be column vector
for i = 1:30
    [jac,f0] = jacobian(func,x);
    if sqrt(dot(f0,f0)/length(x)) < tol
        root = x; return
    end
    dx = jac\(-f0);
    x = x + dx;
    if sqrt(dot(dx,dx)/length(x)) < tol*max(abs(x),1.0)
        root = x; return
    end
end
error('Too many iterations')
function [jac,f0] = jacobian(func,x)
% Returns the Jacobian matrix and f(x).
h = 1.0e-4;
n = length(x);
jac = zeros(n);
f0 = feval(func,x);
for i =1:n
    temp = x(i);
    x(i) = temp + h;
    f1 = feval(func,x);
    x(i) = temp;
    jac(:,i) = (f1 - f0)/h;
end

and printSol :

function printSol(xSol,ySol,freq)

% Prints xSol and ySoln arrays in tabular format.
% USAGE: printSol(xSol,ySol,freq)
% freq = printout frequency (prints every freq-th
%        line of xSol and ySol).

[m,n] = size(ySol);
if freq == 0;freq = m; end
head = '     x';
for i = 1:n
    head = strcat(head,'     y',num2str(i));
end
fprintf(head); fprintf('\n')
for i = 1:freq:m
    fprintf('%14.4e',xSol(i),ySol(i,:)); fprintf('\n')
end
if i  ~= m; fprintf('%14.4e',xSol(m),ySol(m,:)); end

Unfortunately, when I run it : I have an error stating that one of the matrices is singular : I assume that there is somewhere a division by 0 that occurs...

Warning: Matrix is singular to working precision. 
> In newtonRaphson2 (line 23)
In shoot4nl (line 14)
 
Error using newtonRaphson2
Too many iterations

Error in shoot4nl (line 14)
u = newtonRaphson2(@residual,u);

Could someone help me with his mathematical expertise in order to solve this problem ?

Thank you in advance,

Best regards.