For best results, I replaced the constant  2  with a third parameter  a3 . We specify the segment  [a,b]  and the step  h  to obtain data  X  and  Y  and then use the  Statistics:-NonlinearFit  command. When reducing the segment  [a,b] , we get better results:
 

Edit.

Download nonlinfit2.mw

Here is another way. The first and second points almost coincide, so they block one another:

restart;
plot([[[28,.6481496576]],[[28, .648149657615473]],[[28, .6512873548]]],style='point',color=["Blue","Orange","Red"], symbol=solidcircle, symbolsize=12, labels = ["k", "y(k)"], legend = ["10-digit precision", "15-digit precision", "Floating-point iteration"] ,legendstyle = [font = ["HELVETICA", 9], location = right]);

Use single quotes for single use (in 2d math only):

'a >= b'

Example:

-%int(x^2,x);

I think this is not a bug, but just such a design in Maple 2015. In subsequent versions of Maple, this command has been improved (the word  parameters  can be omitted). For example, in Maple 2018.2, the following code works correctly:

y := [1, 3, 8];
val := r->y[r]:
Explore(val(r), r=1..3);

Download diff.mw

In  display(L2,L22,L3,L1)  the previously recorded polygon  L2  closes the later recorded polygon  L22, but L1  and  L3  always lie higher (in Maple 2018.2). See the workaround below:

restart;
with(plots):
L1 := textplot([2, 2, "Polygon"], color = white, font=[times,bold,16]):
L2 := plottools:-polygon([[0, 0], [3, 4], [3, 1]], color = red):
L22 := plottools:-polygon([[0, 0], [0.5, 2], [1,0]], color = green):
L3 := contourplot(x^2 + y^2, x = 1 .. 1.5, y = 4/3*x..2):

display(L2,L22,L3,L1); 

The  coords=polar  option doesn't seem to work (in Maple 2018.2), so I used Cartesian coordinates.
3 regions are plotted:

A:=plots:-inequal({sqrt(x^2+y^2)<2+2*x/sqrt(x^2+y^2),sqrt(x^2+y^2)>3},x=-3.3..4.3,y=-4.3..4.3, color=green):
B:=plots:-inequal({sqrt(x^2+y^2)>2+2*x/sqrt(x^2+y^2),sqrt(x^2+y^2)<3},x=-3.3..4.3,y=-4.3..4.3, color=blue):
C:=plots:-inequal({sqrt(x^2+y^2)<2+2*x/sqrt(x^2+y^2),sqrt(x^2+y^2)<3},x=-3.3..4.3,y=-4.3..4.3, color=red):
plots:-display(<A | B | C>, scaling=constrained);

We can get a general solution to your problem for an integer  n , if we first solve without initial conditions, and then impose these conditions and solve the corresponding system. We see that there is an infinite family of solutions  y(x)=C*sin(n*Pi*x)  (С is an any constant)  only if  A = 0. There are no any solutions if  A <> 0 .

Download ode.mw

It is easy to achieve good visibility by simple means. I changed the style of the surfaces, removing the lines, each plane made in different colors and a few more minor changes. The solution itself is depicted as a bold red dot:

restart; with(plots):
sys := [p+x+.6*y-15, p+.3*x+.2*y-10, p+.5*x+y-14]:
sol:=solve(sys, [x, y, p])[];
A:=implicitplot3d(sys, x = 0 .. 10, y = 0 .. 10, p = 0 .. 10, style=surface, color=["LightBlue","LightGreen","Yellow"]):
B:=pointplot3d(eval([x,y,p],sol), color=red, symbol=solidsphere, symbolsize=15):
display(A,B, axes=normal, orientation=[-20,80], lightmodel=light4);

Should be:

restart;

M:=Matrix(10, 10, [[1, 0, 0, 0, 1/2, 0, 0, 0, 0, 0], [0, 1/2, 0, 0, 0, 0, 0, 1/3, 0, 0], [0, 0, 1/2, 0, 0, 0, 0, 0, 1/3, 0], [0, 0, 0, 1/3, 0, 0, 0, 0, 0, 0], [1/2, 0, 0, 0, 1/3, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 1/3, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 1/4, 0, 0, 0], [0, 1/3, 0, 0, 0, 0, 0, 1/4, 0, 0], [0, 0, 1/3, 0, 0, 0, 0, 0, 1/4, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 1/4]]);

B:=Matrix(10, 5, [[0, 0, 1/3, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 1/4, 0, 0, 0], [0, 0, 1/4, 0, 0], [0, 0, 0, 1/4, 0], [1/2, 1/2, 1, 0, 0], [1, 1/2, 1/2, 1, 0], [0, 1, 1/2, 1/2, 1], [0, 0, 1, 1/2, 1/2]]);

B^%T.M^(-1).B;

Download Collocation.mw

If you work in the real domain then replace  (y^3)^(2/3)  with  surd((y^3)^2, 3) (see help on the surd command for details). Then we see that the result is true for  y<0  as well:

expr:=-1/6*(y^6-6*y^3*ln(x)+9*ln(x)^2)*y^2/surd((y^3)^2,3);
simplify(expr) assuming y>0;
simplify(expr) assuming y<0;

This can be done in many ways. Here are 2 ones:

Download RealSol.mw

We can easily get the answer in a closed form if we use  rsolve  command. To remove the sign of the sum, it is necessary to split this sum into 3 terms:

Edit.

Download seq1.mw