Least_square method for n=p equations/unknown resolutio: Bounds/constraints not respected on proposed solutions

Question

As a Phd Student, working on non-linear problem i try to find solutions, with physical meanings over the following system: Input variables are A, B, Cf, Jp, ks, Pp, R,T. Output variables are : Jw,Js,Cm,Cp,Pf

A = 0.0027                                      #m/s.bar-1
B = 1.19e-08                                    #m/s
R = 0.0821                                      #L.bar/mol.K
T = 25+273                                      #K
Qp = 140                                        #L/h
Jp = (Qp*(1e-3/3600))/(S_membrane)              #m/s
S_membrane= 14                                  #m²
Pp = 0.45                                       #bar
Cf = 7                                          #kg/m3 

def equations(vars):
    Js,Jw,Cp,Cm,Pf = vars  
    return(Jw-A*(Pf -Pp -2*R*T*(Cm*17.11-Cp*17.11)),
            Js- B*(Cm-Cp),
            Cp - (Js*1000)/Jw,
            (Cm-Cp)/(Cf-Cp)-exp(Jw/ks),
            Jp - Js/1000 - Jw)

I expect to solve the problem, considering an least_square approach, initializing with a vector x0 and bounding variables between Xmin,Xmax values as presented below:

x0 = np.array([1e-8,        #Js    
               1e-6,        #Jw
               0.001,       #Cp
               8,           #Cm
               13])         #Pf 

Xmin = np.array([1e-9,      #Js    
               1e-7,        #Jw
               0.000001,    #Cp
               1.01*Cf,     #Cm
               10])         #Pf 

Xmax = np.array([1e-7,      #Js
                1e-4,       #Jw
                0.1,        #Cp
                500,        #Cm
                18])        #Pf

And then, calling the least_squares method:

res = least_squares(equations,x0,method='trf',xtol=1e-13,ftol=1e-13, max_nfev=5000,verbose =1)

print(res.active_mask)
print(res.message)
print(res.success)

It results:

The maximum number of function evaluations is exceeded.
Function evaluations 5000, initial cost 2.1240e+02, final cost 3.3402e-11, first-order optimality 4.10e-01.
[0. 0. 0. 0. 0.]
The maximum number of function evaluations is exceeded.
False
Cf : 7 Cm : 6.964735678220357 Cp : 6.8969269139880875 Jw : -5.395451652744672e-06 Js : -3.721203572887036e-08 Jp : 2.7777777777777775e-06 dP : 57.21883313421709

I dont understand why my system seems to diverge. Also, i printed the value of Jw and Js between iterations that has to be between J,min and J,max and stay positive which is obviously taking negative values. I also tried different "method" such as 'lm' and 'dogbox', output as similar.

Do you have any idea ? Thank you in advance for your time,