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

The output how OP wants:

P:=combinat:-partition(5):
for p in P[1..-2] do
lprint(sort(p,`>`)[]);
od:

  1, 1, 1, 1, 1
  2, 1, 1, 1
  2, 2, 1
  3, 1, 1
  3, 2
  4, 1
 

Edit.

Your system with the specified conditions has no solutions. Since the system does not contain derivatives with respect to  t , we can easily solve it with  dsolve  command, where arbitrary constants  _C1,...,_C4  are some functions of  t . But if we impose the condition you specified on the function  f4 , then as a result we get that the trigonometric function  cos(2*Pi*x)  is expressed as a linear combination exponential functions with real exponents, which is impossible.

Download system.mw

I guess these two pictures (a) and (b)  are probably made in Matlab? Of course, you can do something similar in Maple. In the example below, the first plot from your file was made in Maple 2018.2:

g11:=plot3d(UU(x,t),x=0..xbas,t=0..tbas,color="LightBlue",title='NumericalSolution', titlefont=[times,18]):
xz:=plot3d([x,1,z],x=0..1,z=0..1,style=line,color=blue,thickness=0,grid=[3,6]):
tz:=plot3d([0,t,z],t=0..1,z=0..1,style=line,color=blue,thickness=0,grid=[3,6]):
xt:=plot3d([x,t,0],x=0..1,t=0..1,style=line,color=blue,thickness=0,grid=[3,3]):
display(g11,xz,tz,xt,lightmodel=none,tickmarks=[3,3,6],labels=[x,t,"UU(x,t)"],labeldirections=[horizontal,horizontal,vertical], labelfont=[times,14],axes=frame, view=[0..1,0..1,0..1]);