@Earl 

The procedure ends because only the circles having the radius >= 1/200 are generated.
Note that the line if N > 500 then return was included just for safety, it could be eliminated.
The maths used in the procedure (and some generalizations) is presented in:
Jeffrey C. Lagarias, Colin L. Mallows and Allan R. Wilks,  Beyond the Descartes Circle Theorem,
The American Mathematical Monthly Vol. 109, No. 4 (2002), 338-361.

@mehdi jafari 

When x is in -10000 .. 10000, the oval region inside -6 .. 6  is too small to be seen or even detected by the implicitplot type procedure which is actually used. Note that if you use the suggestion in my answer then the feasible region can be obtained symbolically and explicit, so the plot will be very accurate.

@Ramakrishnan 

eq:=1/(x*y^(2/3))*8.620689655172415*10^(-16)*(-3.11*10^23*x^2*y^(7/6)-3.92*10^19*y^(25/6)+2.14545039999999*10^29*(0.0108*exp(-45.07/y)+exp(-19.98/y^(1/3)-0.00935317203476387*y^2)))/(x+0.015*y^(1.2)):
plots:-inequal(eq>0, x=-10..10, y=0..100, color=yellow);

@mehdi jafari 

You should know better (being your inequality) but usually an inequality has real sub-expressions and e.g. it has y^(7/6).
By the way, usually the domain must by specified (at least when it's not standard).

@ernilesh80 

Once you have a point where the Hessian is negatively defined, you can find a neighbourhood where it is still so.
You can use the minors of the Hessian and evalr.

An example of using evalr:

f:=(x,y) -> sin(x^2-y)*sqrt(x+y)+1.8:
f(1.,2.);
                          0.342529501
evalr(f( INTERVAL(0.9 .. 1.1), INTERVAL(1.9 ... 2.1) ));
              INTERVAL(0.081205985 .. 0.734869569)

So f is > 0 not only at the point (1,2)  but also for x in  0.9 .. 1.1,  y in 1.9 .. 2.1

@ernilesh80 

In my answer I have shown that your function is not concave (for your first set of parameters) and I suggested a way to show ("probabilistically") the concavity for other sets of parameters.

If your function is not concave and want to prove that still a local max is a global max then there are other mathematical notions such as pseudoconcavity. But this is already not Maple related, it is convex analysis.

@ernilesh80 

With DirectSearch or other global optimization package you cannot check the concavity, you simply find the maximum; I understand that this is your goal, and so the concavity is not necessary.
 

I totally agree about 1D/2D and worksheet/document.

I have not installed 2017.1 yet but I noticed that even in 2017.0 the help in
the command line interface (cmaple.exe) has some problems.
In your example with ?diff  some input lines are not displayed
(for example Diff(sin(x), x$n)).

So, probably  this bug propagates to 2017.1 in Standard interface too.
 

@gaurav_rs 

I think that it's too much trouble for such fine-tuning in Maple.
There are much better approaches e.g. TikZ in LaTeX.

About the converse: it is theoretically possible to write a CAS in LaTeX (which is after all a programming language);
of course nobody is going to do this.

Sorry, I do not have this information. I have used it only occasionally, just instead of soving by hand.

BTW, it would be interesting a speed comparison versus Sage.

@marekszpak 

OK, but the algorithm is not very complicated and the code is clean, see e.g.

showstat(LinearAlgebra:-HermiteForm);

I'd suggest to test it before trying to program a better one.

@maple2015 

You must change the BC conditions. As stated the problem cannot be solved because w(x)=0 satisfies the ODE and the BC ; so, for psi you only have the condition psi'(K/2) = psi'(-K/2) = 0, which of course is not enough (the infimum of phi is obviously 0).

@maple2015 

The second objection remains. You cannot speak about "the" integration constant.
For example, the integral of 1/(x+1) can be written as ln(x+1) + C or  ln(2*x+2) + C  or ...
Also, what is this condition when w(x)=sin(sin(x))?

Now about the new formulation. In this case the problem seems to be without solution: taking y=0 the ODE for w(x) has order 2.