You are wrong writting "The BoxPlot command only accepts one list"
See the corresponding help page.

QBoxPlot_mmcdara.mw

@NIMA112 

I won't promote any commercial product: you can easily find some.
Here is an example (likely among many) of a free software FreeFEM

Don't you think it would be much simpler to use a product, commercial or otherwise, capable of solving elliptic equations numerically, rather than thinking about writing FD or FE code in Maple?
Personally, if I wanted to be efficient, I'd choose the first option.

Writting a FD code to solve a Laplace/Poisson equation over a rectangle with Dirichlet bcs on a regular mesh is quite simple.
If the domain is not rectangular it's better fo use FE, maybe on an unstructured mesh... and this becomes quite complex to code even if the methods are well known.
If your PDE contains non linearities or discontinuous parameters a special care to the approximation of the numerical fluxes must be accounted for. here again well known methods do exist, but coding them adds another level of complexity.
I advise you to think twice to the complexity of what you ask. 

Last point: the problem you refer to in your attached file is not well posed, you can't write simultanously  limit(u[2]*(x, y), x = infinity) = 0 and u[1](0, y) = u[2](0, y) : is your domain bounded or not, is  u[1](0, y) = u[2](0, y) a (weak) bc or an internal condition?
If the domain is bounded what is it?

If you want to know which simplification is responsible of the results you get, jdo this

Download which_one.mw

EDITED:
As several simplification methods lead to the same result, you could be interested in knowing those that give the same. Here is a way

a := full_simplify(T[2]):
b := convert(rhs~(a), set):
t := map(u -> lhs~(a[[ListTools:-SearchAll(u, rhs~(a))]]) = u, b):
print~(t):
                                                           (1/2)
        1  /          (1/2)\ /             1              \     
[s17] = -- \50 + 10 41     / |----------------------------|     
        10                   |/      (1/2)\ /       (1/2)\|     
                             \\7 + 41     / \13 + 41     //     

                                   (1/2)              
                             5 + 41                   
         [s3] = --------------------------------------
                             (1/2)               (1/2)
                /      (1/2)\      /       (1/2)\     
                \7 + 41     /      \13 + 41     /     

                                                       (1/2)
           /      (1/2)\ /             1              \     
    [s4] = \5 + 41     / |----------------------------|     
                         |/      (1/2)\ /       (1/2)\|     
                         \\7 + 41     / \13 + 41     //     

                                       1  (1/2)
                [s9, s8, s7, s6, s5] = - 2     
                                       2       

        [s16, s15, s14, s13, s12, s11, s10, s2, s1] = 

                            /      (1/2)\               
                         10 \5 + 41     /               
          ----------------------------------------------
                           (1/2)                   (1/2)
          /          (1/2)\      /           (1/2)\     
          \70 + 10 41     /      \130 + 10 41     /     

@Axel Vogt 

It's not really frustration, because I don't think my comment could have caused a scientific revolution, and I didn't spend more than 10 minutes writing it. It's more a question of exposing a practice that seems to me to be more and more common.

@MaPal93 

It is likely that I had some script en betwwen the definition of pts and the NLPSolve command.
Whatever I executed again the file  display_formula_mmcdara.mw  and I got the same results as yours: display_formula_mmcdara_2.mw

Sorry for the inconvenience.

By the way I sent you a new reply here

Here corrected:

 

C2 := k>=Q1, Q1>=rho*k   # k>=Q1>=rho*k; IS NOT an inequality, write C2 as a sequence of inequalities
C3 := I2 >= I1:          # Strict inequalities ARE NOT supported

constraints := eval({C1, C2, C3}, DATA);  # inequalities MUST BE evaluated too

# Use Optimization:-NLPSolve to account correctly for the constraints

Optimization:-NLPSolve(TRC(I1, I2), constraints, I1 = 0 .. 5, I2 = 0 .. 7, assume = nonnegative)
     [7.515*10^9, [I1 = 1.500, I2 = 2.500, Q1 = 1.200*10^8, Q1a = .4600, Q1b = 0., Q2 = 5.000*10^8]]

I believe that a simple analysis is enough to get this result (look to my comments to your previous question)

@MaPal93 

opt[1]  is the value of the minimum of the objective function J.
Given the definition of J, opt[1]  is sum of squares of the residual and is exactly what

Statistics:-NonlinearFit(K(Gamma), convert(pts, Matrix), Gamma, output=residualsumofsquares)

 would return if it was capable (Maple 2015) the find the best fir (which would be opt[2])
As a rule I prefer using Optimization:-NLPSolve instead of Statistics:-NonlinearFit: the former is more vesatile (you can account for ranges, initial points and your own expression of the objective function) and less often fail than the latter. 


You wrote "Can you please help me? (I didn't find something like this)"
I believe the attached file could help you. It contains a procedure named NonLinearLeastSquares which wraps a  Optimization:-NLPSolve in order to solve non linear least squares problem.
The way to use it mimics Statistics:-NonlinearFit's

Download NonLinearLeastSquares.mw

Ice on the cake: NonLinearLeastSquares_and_Kriging.mw completes the previous file by adding a simple implementation of Kriging.
The "approximation function" is an exact interpolator which passes through points pts.  Its approximation error out of these points is less than 1e-7.
Keep in mind this Kriging part of the file is purely illustrative and serves only to demonstrate the power of the Kriging approach (note that the interpolator depends on a single parameter named psi whose meaning is a correlation length, that is how far from a point x[i] the information it contains (y[i}) "propagates" to affect the valur y[k] at any other point x[k]).

@MaPal93 


 

1. If I am not satisfied with the linear fit...
   There exist very versatile interpolation models that could be used instead of more classical nonlinear least squares.
   I thing more specificaly to kriging models.
   The CurveFitting package proposes such a model since Maple 2019, proposes an implementation, but having had a look at it, it's not really suited to modeling the output of a computer code.
   Take a look at this approach and let me know if you're interested.
    
2. I guess that using option grid=[N, N] with N larger than the default value (around 30 if I'm not mistaken) should generate more points on the level curve, but be aware that the computation time probably varies like N^2.
    
3. The explanation is provided in the last part of file slices_of_three_surfaces_mmcdara.mw

@MaPal93 
 

1. I didn't check the value, but yes, the Lennard-Jones potential has an asymptote.
2. The green curve is indeed the graph of the potential.
   And yes, its expression is

   eval(K(Gamma), opt[2])

3. Yes, ar least iw we consider que domain 0 < Gamma < 10, -1 < rho < 1.

"Why don't you embed the L-J snippet in the worksheet and post your comment as an answer I can upvote/give best answer?".
Because I'm fed up with this voting and ranking system.
I'm not participating to gain notoriety but to provide help and share my limited experience. That is why I decided to publish only comments and satisfy myself with a simple positive feedback. 


Error, (in CurveFitting:-Spline) data points not in recognizable format
I use Maple 2015. 
Here is the detail of the operators in pc0:

# with Maple 2015 op(pc0) returns something like this:

CURVES(
    [[x[1], y[1], [[x[2], y[2]],
    [[x[3], y[3], [[x[4], y[4]],
    ...
    [[x[n], y[n], [[x[n+1], y[n+1]],
    COLOUR(RGB, .47058824, 0., 0.54901961e-1)
),
AXESLABELS(GAMMA, rho)

# so op(1, pc0) returns 

CURVES(
    [[x[1], y[1], [[x[2], y[2]],
    [[x[3], y[3], [[x[4], y[4]],
    ...
    [[x[n], y[n], [[x[n+1], y[n+1]],
    COLOUR(RGB, .47058824, 0., 0.54901961e-1)
)

# and [op(1..-2, op(1, pc0))] returns 

[
    [[x[1], y[1], [[x[2], y[2]],
    [[x[3], y[3], [[x[4], y[4]],
    ...
    [[x[n], y[n], [[x[n+1], y[n+1]]
]

It is higly likely that the structure of pc0 that your newer version builds is not the same.
Can you execute the command op(pc0) in yourworksheet and send it to me?

Here is an alternative to fsolve for the "Surface 1" case only.

Explanation:

1. plot the level curve  H(Gamma, rho) = 0 where H(Gamma, rho ) = D[2](f)(Gamma, rho),
2. extract the points (Gamma[i], rho[i]) (remove duplicate abscissas),
3. build the spline approximation of the curve passing through points  (Gamma[i], rho[i]): this curve is a piecewise representation of the function Z : Gamma --> Z(Gamma) such that H(Gamma, Z(Gamma)) = 0.
    

Download Surface1_fsolve_workaround_mmcdara.mw

n case an approximation of the function Z(Gamma) would be enough for you, you can observe that the level curve   H(Gamma, rho) = 0 looks like a Lennard-Jones potential.
Then you can try to fit it this way:

K := Gamma -> a+b*((c/Gamma)^d-(e/Gamma)^f):
J := add((map2(op, 2, pts) -~ K~(map2(op, 1, pts)))^~2):
opt := evalf[6]( Optimization:-NLPSolve(J, assume=nonnegative) );

display(
  plot(pts, style=point, symbol=circle, color=blue), 
  plot(eval(K(Gamma), opt[2]), Gamma=1..10, color=green, view=0.5..1)
);
opt := [0.164888e-2, [a = .747049, b = .980991, c = 1.18492, d = 1.89732, e = 1.14613, f = .985049]]

Here is the result.

You already wrote a quite impressive Maple code (or maybe did you copy it from some source?)

If you indeed write the cote which implements the N-soft set method, I suppose you know what this method is and how to use it (personally, I don't code a method I don't know about or I don't don't understand the results ir provides.).
If you know what the N-soft sets method is you shouldknow what results it provides andhow to analyze them, don't you?
 

Maybe you should use a method you know how to analyze the results.

In any case, your question isn't strictly a Maple question (you've already written (?) the code) but a Statistics question.
 

@MaPal93 

I try to find something explaining the results you get with ncal3 but I failed.
Probably because I didn't know what I was seeking...?

By the way, your code stops at the "discontinuity" with Maple 2015.

@Andiguys 

You're welcome

@FDS 

M[2..-1] -~ M[1..-2]

But if you prefer a procedure which returns a vector:

tst := proc (M) 
   local N, i, m; 
   N := numelems(M):
   m := Vector(N-1):
  for i from 2 to N do 
     m[i-1] := M[i]-M[i-1] 
   end do; 
   return m
 end proc:

tst(M)