MATLAB resolution of a 6-th order non-linear differential equation [stiffness]

Question

I am currently trying to solve a non-linear differential equation of order 6 for a function F defined on :

Or :

With these following boundary conditions :

And an additional integral condition that is written :

I am using MATLAB in order to do so with the "BVP5C" solver. I wrote my integral constraint as a new ODE that I implement in the solver.

Here is the script :

%% Solving the problem of elastic droplet's spreading in an elastic regime in the [0;eta0] interval 

% Solving method using bvp4c solver
% Boundary conditions : F'(0) = F'''(0)= F'''''(0) = 0 and F(eta0) = F'(eta0) = F''(eta0) = 0
% Integral constraint on the volume

eta0 = 4;
yinit = [1; 0; 0.5079; 0; 0.1194; 0; 1]; % Initial conditions for the bvpinit solver.
nbiter = 10; % Number of iterations.

% We choose to iterate in a log scale the values of eta
% We integrate by using the bvp4c solver for the corresponding ODE's

eta0tab = logspace(0,log10(eta0),nbiter);


for k=1:nbiter
    eta = linspace(0,eta0tab(k),100);
    solinit = bvpinit(eta,yinit);
    sol = bvp5c(@fun_ode, @bcfun, solinit);
    eta = sol.x;
    F = sol.y(1,:);
    yinit = sol.y(:,1);
end

% We plot the result [F; F'] in a new figure

figure()
plot(eta,-F,'-','LineWidth',2);
grid on
hold on
plot(eta,sol.y(2,:),'LineWidth',2)
legend('F',"F'",'Location','northeast')
xlabel('\eta','FontSize',13);
ylabel("F(\eta),F'(\eta)",'fontsize',13)

% We first construct a vector of 1-st order ODE's
% The last one is the constraint term

function dF = fun_ode(eta,F)
    
dF=[F(2:6);
    ((1/(9))*(5*F(1)-4*eta*F(2))-3*F(1)^2*F(2)*F(6))/F(1)^3; F(1)];
end

% We built residuals of the method with y(b)-a = 0.

function res = bcfun(ya,yb)
   
res =[ya([2 4 6]); yb(1:3);yb(7)-1];

end

Unfortunately, as I am changing the value of the upper boundary of the interval I keep having different shape of the solution.

My question is :

Is it an error of my coding ? Or is my problem not stiff at all with these boundary conditions ?

For instance, here is the result for : eta0 = 5 :

Then for eta0 = 15 :

Thank you in advance for your help.

N.B : Since there is an integral constraint, I had to remove a boundary condition :