error bvp4c a single jacobian encountered

Hi, I'm trying to plot the P value (corresponding to a peak in a P-umax plot) in function of "s", a variable linspaced between 0.1 and 2.

I've tried to write something like this but I don't know what to put in the Pcr findpeak function and I'm obtaining an error for the bvp4c.

I've tried to write this code:

clc;
clear;
clear;
%------------------------------------------------------------------------
ni=0.3;
L=50;
R=20;
E=210000;
E1=E/(1-ni^2);
%-----------------------------------------------------------------------
s= linspace(0.1,2,10);
I = zeros(size(s));
k = zeros(size(s));
Pcr = zeros(size(s));
%%
P=linspace(7800,3300000,10);
y_max_values=  zeros(size(P));
for i=1:length(s)
    for j=1:length(P)
     f = @(z,u) [u(2); u(3); u(4); (-1/(E1*I(i)))*((P(j)*u(3))+(k(i)*u(1)))];
     bc = @(ua, ub) [ua(1); ua(2)-0.00000001; ub(1); ub(2)+0.00000001];
    zmesh = linspace(0, L, 100000);
    iguess= @(z) [0; 0; 0; 0];
    solinit = bvpinit(zmesh, iguess);
    Sol= bvp4c(f, bc, solinit);
    
    x = Sol.x;
    y = Sol.y;
    idx = (x > 50) & (x < 180);
    x_new = x(idx);
    y_new = y(1, idx);
    
     y_max_values(j) = max(y_new);
    Pcr=findpeak (?) %what I should put here?
    end
    
end
%%
figure
hold on
grid on
plot(s, Pcr)
xlabel ('s [mm]')
ylabel ('Pcr [N]')

I'm obtaining this error:

Error using bvp4c (line 251)
Unable to solve the collocation equations -- a singular Jacobian encountered.

Error in plot_s_Pcr (line 33)
    Sol= bvp4c(f, bc, solinit);

I've written a code that plots P-umax for s=1 but I need it to be repeated for every s value

clc;
clear;
close all;
R=20;
L=3*75;
s=1;
ni=0.3;
E=210000;
w=2;
K=0.5;
%
E1=E/(1-(ni^2));
I=(2*pi*R*(s^3))/12;
k=2*pi*((E*s)/R);
P = linspace(400000,9000000);
y_max_values = zeros(size(P));
for i=1:length(P)
    f = @(z,u) [u(2); u(3); u(4); (-1/(E1*I))*((P(i)*u(3))+(k*u(1)))];
    bc = @(ua, ub) [ua(1); ua(2)-0.00000001; ub(1); ub(2)+0.00000001];
    zmesh = linspace(0, L, 100000);
    iguess= @(z) [0; 0; 0; 0];
    solinit = bvpinit(zmesh, iguess);
    Sol= bvp4c(f, bc, solinit);
    
    x = Sol.x;
    y = Sol.y;
    idx = (x > 50) & (x < 180);
    x_new = x(idx);
    y_new = y(1, idx);
    
    y_max_values(i) = max(y_new);
end
figure
hold on
grid on
plot(y_max_values, P)
legend ('cilindro 1 w0')
xlabel ('u_max [mm]')
ylabel ('P [N]')

So I need to find in the plot, the P value corresponding to the first peak of the following P-umax plot, where the umax value moves from zero in the first time. (8x10^5 in this case), and repeat this procedure for every "s".

The P value corresponding to peak I'm trying to find has to be assigned to "Pcr"

Could someone help me out ?

Is the the code structure right for the aim I need to?