@ecterrab 

I find myself in your comment I agree with completely.

@vv 

When I receive Answers or Replies to my questions I read all of them all just as carefully whatever their status.

And, okay, in the Mapleprimes logic I sometimes though "This Comment deserved to be an Answer, for I would have rewarded you for it and thus raised your value/credibility for the future".
I think that a lot of people are happy to be rewarded this way, and yes, I was too until a while ago. Until I sensed a kind of falseness behind this reward system and simply decided not to enter the game any more.

As I replied to @acer  OP's opinion is more valuable to me than the number of badges I collect.
@acer  raised the issue of Reply appearing before Answer, thus rejecting at the end of the list the potentially more elaborated contributions.
This can be changed programmatically or, I can send only answers with a header telling "Do not vote for this answer nor select it as the best". 

@acer 

for having changed  "MapleSim" to "MaplePrimes", sort fir the inconvenience.

Yes, I have indeed noticed that a Reply (or a Comment, because when you edit an Answer, there is no "Convert to Reply" option, but a "Convert to Comment" option) appears first of an Answer and I always found this quite curious.
It seems to me the opposite would have been more natural.

To be clear: I'm not doing comments instead of replies to be at the top of the list the OP is about to read, but because I don't care about being rewarded for the work I've done.
I consider an OP answering "very useful, thank you" is much more gratifying than a few points awarded by I don't know who.

@Carl Love 

• Can I confirm that the maxima are global (NLPSolve tends to return local optima)?
  No.
  One often says that starting from many different initial points and ending at the same point gives a strong evidence that this point is indeed a global optimum. But even in this case one cannot be one hundred percent sure of that. Idealy one should use a global optimizer to be conclusive.
  
  More than this I am not even sure iterationlimit is large enough to consider the results are "converged".
  Last point: the solution I get for n=20 doesn't look like the solution the OP presents. It would br interesting to know what sum of radii this latter corresponds to.
   
• That's quite impressive work!
  Thanks, but I'm not concinced of that either.
  I just modeled the problem in a simple way and used Maple do find (one of) its solution(s). Nothing impressive behind this.
  
  In case you would be interested in packing problems here is a erference site packomania which gathers a lot of results concerning optimal packings.

@spalinowy

Thank you for selecting my answer as the best one (there was no other anyway :-) ).
However, I'm going to convert it into a comment for the very simple reason that I no longer wish to participate in the Mapleprimes rewards system.
Doing this will make me lose a few points, but the only thing I value is that my contribution has helped you, which seems to be the case.
All the best

WHOEVER YOU ARE, PLEASE STOP TURNING THIS COMMENT INTO A REPLY.
THANKS IN ADVANCE



Check this

Download Disks_within_unit_square.mw

 

WHOEVER YOU ARE, PLEASE STOP TURNING THIS COMMENT INTO AN ANSWER.
THANKS IN ADVANCE


If you want to get a NUMERIC solution, please look to this file.

There are several points that it is up to you to fix:

1. There are two functions

    {z[F0](t), z[R0](t)}

   that are not defined. 
   If you want to solve nimerically the sustme you must give them an analytic expression.
   I made an arbitrary choice to go further.

2. You missed 4 initial conditions.
   I made an arbitrary choice to go further.

3. Either you give the parameters a numeric valur OR you use(this is what I did) the "'parameters' option
   I made an arbitrary choice for their values to gget a fully instanciated solution

Download How_to_get_a_solution.mw

@Aung 

Thanks for your feedback.
Feel free to ask for more help on this problem.

By the way: the Statistics:-NonlinearFit solver uses Optimization:-NLPSolve's and this latter supports constraints (for instance p3 > p2, p1 > 0, ...) which makes it, in my opinion, more versatile and powerful than the former, to solve non linear least square problems.
The main differences between Statistics:-NonlinearFit and Optimization:-NLPSolve are:

• the former builds the function to be minimized (the summ of squared residuals) while you must define it explifitly with the latter;
   
• the former provides more outputs than the latter... but it is very easy to reconstruct them as soon as Optimization:-NLPSolve has privided the location of the minimizer.

@acer 

@acer 

Thanks for the information about loaded packages.

@acer 

Before the 2nd Tabulate you shoulf change opts:

opts := fillcolor = (proc (T, i, jj) options operator, arrow; `if`(jj = 1, "LightBlue", "PeachPuff") end proc), widthmode = pixels, exterior = none:

By the way, the OP loads Student:-Statistics and Statistics: can this generate conflicts?
For instance, if one calls Mean, which one is used?

I'm admittedly a bit naive, but I had hoped to get to have some feedback my answer to your previous question.
Whatever I see you understood how to build an unweighted sample from a weighted one, so at least you read my reply...

By the way (1) , the formula you use to compute EcartType is uncorrect, the denominator should be add(Weights)-1

By the way (2), still another Maple Bug:
StandardDeviation(data, 'weights'=Weights) and Variance(data, 'weights'=Weights) return wrong results.

don't you ask your question here MMA ?

@MaPal93 

Then what do you mean by "the interpolator depends on a single parameter named psi"?
The command to construct the kringing emulator is

krig := KG_Matern52(data[[seq(1..59, 1)], ..], 5):

where 5 is THE SINGLE parameter psi the solution depends on (compare to the 6 parameters of the Lennard-Jones model).

Feel free to change this value, for instance use 20 instead of 5: you will get a curve which still passes through the points (its is an interpolating curve) but which strongly oscillates between these points.
The value of psi can:

• either be set to a convenient value trom a visual observation of the resulting interpolation curve (what I did),
• or searched by solving a specific minimization problem (which I didn't implemented here).

The expression is indeed complex and is never displayed. A kriging emulator is always used as a black-box function: you give it the value of the input(s) and it returns the corresponding value of the output (and possibly some other useful informations).
To draw a parallel, Its usage is exactly the same as sol in

sol := dsolve(ode_system, numeric)

Kriging is an extremely powerful and versatile method which is widely used in M&S (Modelling and Simulation) to build analytic models (aka metamodels) of outputs of a costly computer code.

For insance your Maple code can be seen as a complex function F : Gamma_1 --> Gamma_2 = F(Gamma_1) whose no closed-form is known.
Observing a collection C = { (Gamma_1[n], F(Gamma_1[n])),  n=1..N } of points is sometimes enough to infer a possible expression for F.
But not always, for instance I'm not sure that a lot of people would have say here "A correct metamodel is a Lennard-Jones type model".

The main problem with classical metamodels (a linear one for instance), is that there is almost no chance that it passes through all the points C contains.
How can you justify this when F is a deterministic function and when there is no observation noise on the outputs Gamma_2?
(The situation is completely different when points come from empirical (aja experimental) observations).

So you have to build an interpolator and not an estimator, which means that if M is the metamodel, you want M(Gamma_1[n]) = F(Gamma_1[n]) for all n=1..N.

There exist an infinity of such estimators (probably the most known are spline interpolators of any odd order), Krining (aka Gaussian Emulator, GE in the sequel) is another type.
The simplest GE depends on a single parameter (I denoted psi) for scalar inputs. For vector inputs psi is then a vector of same size. Beyond this simple GE there exist more complex variant dependong on parameters whose natures differ from psi's.
To be simple psi characterizes the smoothness of the interpolator. If you believe the interpolator must be very regular, or smooth, you must use a "small" value for psi, a "large"one if you believe the interpolator can have large deviations out of the points in C.

The fact that the smoothness of the GE is almost the only parameter to fit makes Kriging one of the most powerful tool to build high dimensional metamodels (in dimension d psi contains d parameters: compare this number to the fit a polynomial of total degree q in dimension d which depends on (q+d)!/q!/d! parameters) 

At the end the GE is always used as a black-box function: you give it the values of th inputs (Gamma_1 in your case) and it returns the output (here Gamma_2).
Despite the fact it has a closed-form expression this latter can be that complex that it is never displayed (that is exactly the same thing for a cubic spline interpolator for instance).


Unfortunately there is no easy to read synthetic reference about gaussian emulators
In coase you go further I can browse the web to try and find a few.

@janhardo 

An approximation is not a simplification.
If you think the opposite, there are two even simpler expressions:

taylor(expr1, x=0, 5);
                            31  3    / 5\
                      2 x - -- x  + O\x /     (basically what S22 gives)
                            3            
taylor(expr1, x=0, 1)
                              O(x)   (4 characters alone)
mtaylor(expr1, x=0, 1)
                               0   (only one single character: you can't do simpler)