@JAMET 

If it so look at this Youtube video.
The first conic the speaker talks about is a parabola (your case).
Follow step by step what he does:

1. use LinearAlgebra:-Eigenvectors to diagonalize the Matrix(2$2, [9, -12, -12, 16]) associated to your quadratic form, let P the matrix of right eigenvectors, use the transformation <u, v> = P^(-1).<x, y> to align the parabola axis with u, (or v see below) ;
    
2. use Student:-Precalculus:-CompleteSquare to get the expression a*(u-p)^2 + b*v + c (or a'*(v-p')^2 + b'*u + c' depending on the order LinearAlgebra:-Eigenvectors returns the eigenvalues);
    
3. rewrite b*v + c into b*(v + c/b)  (or b'*(u + c'/b')) to get the translation aling u (v);
    
4. set U = u-p and V=v-c/b (or V = v-p' and U=u-c'/b' ) to get the reduced form a*U^2+b*V=0 (or a'*V^2+b'*U=0)
    

@dharr 

If I have some time I will try to design a more subtle algorithm.
 

@schachar 

You wrote " Please use the convolution of field (X, Y) by a 2D  gaussian for the data", but you did not provide any field.
You provided two series X and Y index by n (n=1..781).
A discrete 2D field F(X, Y) is a set of triples (X_i, Y_i, F(X_i, Y_i)). There is no F(X_i, Y_i) here.

@vv 

The epsilon was intended to be a workaround to the conditions r[i] > 0.
My intention was to write instead r[i] >= epsilon, that I changed to r[i]+epsilon >= 0 by carelessness.

Thanks for pointing out this stupid mistake.

@minhthien2016 

I got 1.905667204, so almost the same thing.

Do you know what algorithm your friend used (the (apparent?) symmetry with respect to the second diagonal is remarkable)?
Are solutions obtained by Asymptote for other values of n still symmetric with respect to a diagonal?
Those I got (see  mmcdara's reply) are symmetric for n from 1 to 9 only.

Meanwhile I edited my reply because the plots were not correct. You will se that the plots are now identical to the one you get with Asymptote.
 

By the way sand15 / mmcdara is one and the same person (provessional/personal login)

@acer 

I didn't check.
My intuition is that there should be equalities.

The situation would be different if you search the best way to pack N disks or same radius R into the unit square. In this case they may exist somewhere a free volume of size > Pi*r^2 where a last disk can be placed in an infinity ways, not even touching the other ones nor the boundary of the box.
But if you allow this last disk to have a radius R' > R to "fill the void", then there exist only one posibility to place it when R' reaches some maximum value.

 

@vv @Carl Love

EDITED
 

Some digging on the n=14 case.
I ran 3 sequences of 1000 optimizations all starting from a different initialpoint randomly chosen (unless if I am extremely unlucky all those initial points respect the constraints).
For these three sequences I kept the packing which gives the highest sum of radii (S).

I got three times the same value S=1.90567, the same sets of radii, and the packings are identical up to rotations by an angle multiple of Pi/2.

I belive these results partly answer @Carl Love's question "Can we confirm that the maxima are global?".
Indeed, I would infer that S=1.90567 is a global maximum.

As the two histograms provided seem to demonstrate, there could exist a large number of local maxima where NLPSolve is stucked.

 

>  restart >  with(plots): with(plottools): with(Statistics): DISKS WITHIN THE UNIT SQUARE The packing shiuld be such that the sum of the radii is maximum >  epsilon := 1e-8: n := 14: J := add(r[i], i=1..n): P := [seq(<x[i], y[i]>, i=1..n)]: cons_1 := seq(seq(add((P[i]-P[j])^~2) >= (r[i]+r[j])^2, j=i+1..n), i=1..n-1): cons_2 := seq(op([P[i][1]-r[i] >= 0, P[i][1]+r[i] <= 1, P[i][2]-r[i] >= 0, P[i][2]+r[i] <= 1]), i=1..n): cons_3 := seq(r[i] >= epsilon, i=1..n): SoR : Sum of the Radii Max_trials : maximum number of trials                      Note: I used a while loop instead of a for loop because NLPSolve sometimes returns an error >  Packing := proc(n, Max_trials)   local SoR, SoR_max, error_found, error_count, rep, StartAt, opt, sor, error_type, Best_packing:   SoR     := Vector(Max_trials):   SoR_max := 0:   error_found := false:   error_count := 0:   rep := 0:   while rep < Max_trials do     StartAt := {                  seq(x[i]=rand(0. .. 1.)(), i=1..n),                  seq(y[i]=rand(0. .. 1.)(), i=1..n),                  seq(r[i]=rand(0. .. 1e-6)(), i=1..n)                };        try       opt := Optimization:-NLPSolve(J, {cons_1, cons_2, cons_3}, maximize, iterationlimit=10^4, initialpoint=StartAt):     catch:       error_found := true     end try:        if `not`(error_found) then       rep := rep+1:       sor := evalf[6](eval(add(r[i], i=1..n), opt[2]));       SoR[rep] := sor:          DocumentTools:-Tabulate(         [sprintf("%d", rep), sprintf("%1.6f", SoR_max)],         'alignment'='left',         'widthmode'='percentage',         'width'=20       ):              if sor > SoR_max then         SoR_max      := sor;         Best_packing := opt[2];       end if:     else       error_found := false:       error_count := error_count+1:       error_type  := lastexception:     end if:   end do:   return  table([             Errors = `if`(error_count = 0, None, [Type = error_type, Occurrences = error_count]),             RadSum = SoR,             Best   = Best_packing           ]) end proc: >  res := Packing(14, 1000): res[Errors]; SoR := res[RadSum]: Best_packing := res[Best]: Best_packing_1 := copy(Best_packing): (1) >  FivePointSummary([SoR]): Histogram([SoR], minbins=100) >  interface(rtablesize=n+2): `<,>`(   `<|>`(['i', 'x[i]', 'y[i]', 'r[i]']),   convert([seq(eval([i, x[i], y[i], r[i]], Best_packing), i=1..n)], Matrix),   `<|>`([`sum`, ` `, ` `, eval(add(r[i], i=1..n), Best_packing)]) ); (2) >  Histogram(   (evalf[6]@rhs)~(select(has, Best_packing, r))   , discrete=true   , size=[1200, 100]   , color=red   , thickness=3   , title="Radii" ); Packing_1 := display(   seq(     disk( eval([x[i], y[i]], Best_packing), eval(r[i], Best_packing), color=blue)     , i=1..n   ),   rectangle([0, 0], [1, 1], color="LightGray")   , title=typeset(Sum(r[i], i=1..n) = evalf[6](eval(add(r[i], i=1..n), Best_packing)))   , size=[500, 500]   , scaling=constrained ): RUN 2 >  res := Packing(14, 1000): res[Errors]: SoR := res[RadSum]: Best_packing := res[Best]: Best_packing_2 := copy(Best_packing): >  FivePointSummary([SoR]): Histogram([SoR], minbins=100): `<,>`(   `<|>`(['i', 'x[i]', 'y[i]', 'r[i]']),   convert([seq(eval([i, x[i], y[i], r[i]], Best_packing), i=1..n)], Matrix),   `<|>`([`sum`, ` `, ` `, eval(add(r[i], i=1..n), Best_packing)]) ): Histogram(   (evalf[6]@rhs)~(select(has, Best_packing, r))   , discrete=true   , size=[1200, 100]   , color=red   , thickness=3   , title="Radii" ); Packing_2 := display(   seq(     disk( eval([x[i], y[i]], Best_packing), eval(r[i], Best_packing), color=blue)     , i=1..n   ),   rectangle([0, 0], [1, 1], color="LightGray")   , title=typeset(Sum(r[i], i=1..n) = eval(add(r[i], i=1..n), Best_packing))   , size=[500, 500]   , scaling=constrained ): RUN 3 >  res := Packing(14, 1000): res[Errors]: SoR := res[RadSum]: Best_packing := res[Best]: Best_packing_3 := copy(Best_packing): >  FivePointSummary([SoR]): Histogram([SoR], minbins=100): `<,>`(   `<|>`(['i', 'x[i]', 'y[i]', 'r[i]']),   convert([seq(eval([i, x[i], y[i], r[i]], Best_packing), i=1..n)], Matrix),   `<|>`([`sum`, ` `, ` `, eval(add(r[i], i=1..n), Best_packing)]) ): Histogram(   (evalf[6]@rhs)~(select(has, Best_packing, r))   , discrete=true   , size=[1200, 100]   , color=red   , thickness=3   , title="Radii" ); Packing_3 := display(   seq(     disk( eval([x[i], y[i]], Best_packing), eval(r[i], Best_packing), color=blue)     , i=1..n   ),   rectangle([0, 0], [1, 1], color="LightGray")   , title=typeset(Sum(r[i], i=1..n) = eval(add(r[i], i=1..n), Best_packing))   , size=[500, 500]   , scaling=constrained ): The three best solutions found are identical after ad hoc rotations >  DocumentTools:-Tabulate([   Packing_1,   translate(rotate(translate(Packing_2, -1/2, -1/2), -Pi/2), 1/2, 1/2),   translate(rotate(translate(Packing_3, -1/2, -1/2), -Pi), 1/2, 1/2) ]) >

Download Disks_within_unit_square_14_edited.mw
 

@mehran rajabi 

Is the expression of f given once and for all or is it just a notional example that you provide?

Do you want Something_like_this.mw ?

Here are A_few_generalizations.mw

Run Maple, open a new worksheet or document and execute this command

help(dsolve)

Read the help page and learn how to solve a differential equation with Maple.
Hint:

# y' = y+x means y'(x) = x+y(x) and is written diff(y(x), x) = x+y(x)
#
# At this point there is no need to use complex stuff like declaring that y is an alias for y(x) or
# using special tricks to us y' instead of diff(y(x), x)

dsolve(diff(y(x), x) = x+y(x))
                   y(x) = -x - 1 + exp(x) _C1

upload your worsheet by using the green arrow in the menu bar

UPDATED

This file contains the solutions of 3 well-posed problems:

1. {f(0) = 0, (D(f))(0) = 0, ((D@@2)(f))(0) = a, theta(0) = 1, (D(theta))(0) = b}

2. {f(0)=0, D(f)(0)=0, f(10)=1, theta(0)=1, theta(10)=0}
   

3. {f(0)=0, D(f)(0)=0, D(f)(10)=1, theta(0) = 1, theta(10) = 0}
   

Make your choice.
And if none of these three problems corerspond to what you want, try to take inspiration of what is done in the attached file to solve your problem of interest.

Several_well-posed-problems_sand15.mw

@acer 

Do you think it's decent to improve your reputation based on work others have done?
You don't need to. You've already earned it.
Personally, I would have mentioned Kitonum's work in a simple reply, not an answer.

@janhardo 

As ChapGPT is meant to provide clear synthesis on a subject, the least one can hope for is that its answers are error free and understandable.
For instance

• Can you explain me how ChatGPT passes from the 4th equality to the 5th one (in point 4)?
• Can you explain me how ChatGPT finds that (point 7) pi³/3! = pi²/6?
• More generally, ChatGPT does not shine for its clarity (Point 6 is quite obscure).

Better sources than CheatGPT exist  to understand how Euler solved the Basel problem.
Wiki (see @vv) is already better but nothing is worth the original  Euler (english translation)
 

@vv 

Thank you for this explanation (even if I only remember the analytical continuation in name, it's all too far away for me).
If I have a bit of time, I'll try to refresh my memory.

@Jesús Guillera 

By the way, did you browse the references I provided you or did you made your own research about convergence conditions of your serie?

To end my contribution: are you sure what you said in your last reply is correct?
 

kernelopts(version);
Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895

F(-2);
                   infinity                 
                    -----                   
                     \              2     n 
                      )     po(1, n)  (-2)  
                     /     -----------------
                    -----          2        
                    n = 1     /1   \         
                           po| -, n |  (n + 1)
                              \2   /         

A truncated serie:

FT := (z, N) -> add(po(1,n)^2/po(1/2,n)^2*z^n/(n+1),n=1..N):
evalf(F(1/2)), evalf(FT(1/2, 100)); 
                    2.254576561, 2.254576561

 evalf(F(-2)) , evalf(FT(-2, 10)), evalf(FT(-2, 50)), evalf(FT(-2, 250)); 
                                          15                75
 -1.076461643, 2002.463620, 2.323658088 10  , 3.777981621 10  

Convergence doesn't seem to hold at z=-2...

Look at this: Download Convergence_criteria.mw