add the option range=-1..1 as the second argument to HeatMap works for me in Maple 2017.

The 'ps' driver (not 'eps) is supposed to be uncolored by default, but didn't seem to be. Using color=none is black and white,  but for me it produces a black background with white text and lines. Using color=gray gave a gray plot line, but no axes or text.

So adding color="Black" in the display command to override default plot colors should be an alternative, but if colors were given in the plot command, then you probably have to edit the plot structure.
 

The general solution can be (partly) obtained as below. Here sqrt(_c[1]) is n, but to get this you need to know that the function is periodic in theta (trying u(r,theta)=u(r,theta+2*n*Pi) doesn't work). To keep the solution finite at the centre of the disk you have to set _C2=0, and again you have to manipulate with Maple further to decide this.

Download PDE2.mw

Using HINT=`*` forces the product form, and then you can implement another BC. Edit: added the 2nd BC "by hand"

Download PDE.mw

As the others have implied, a lot depends on what you consider simple. For the last one I like

Student[Precalculus]:-CompleteSquare(k^2+2*k+9); which gives (k+1)^2+8

and for the rational functions I like continued fractions:

convert((k-2)*(k^2+5)*(k^3-k^2+7*k+8)/(6*k*(k^2-3*k+8)),confrac,k);

The second sqrt argument k^4-10*k^3+37*k^2-60*k+180 can be simplified to (k-2)^2*(k-3)^2+144. The first one is harder.

Download sqrts.mw

Your answer to @Preben Alsholm  suggests you are using a table, not an array. Then the following works:

Download Occurrences.mw

Under help for list (in Maple 2017.3), we find:
"Lists and sets can be nested, in which case selection can be done in one of two ways:  S[i][j]...[n] or S[i,j,...,n]."

This seems to allow the interpretation of m[[1,2],2] as m[[1,2]][2] though it is very confusing.

The use of f(x) and df(x) is confusing in this context. And the print is redundant. Just use:

restart;
f:=x^2-3;
df:=diff(f,x);
x:=1.;  #use floating point here
for i from 2 to 5 by 1 do
    x:=x-(f/df);
    end do;

Your second approach works if you use solve:

Download linear.mw

If I paste your code into an execution group and hit enter, I see the cursor just before R on the line starting: locals R (It's quite hard to see). This will be near the error, which is that "locals" should be local. Continuing and fixing some missing semicolons gives a syntax error free code:

 # ctrl + del to delete a Maple cell
 # golden search implementation
 # chapra 7th ed
 golden_search := proc(f, xl, xu, es100, maxiter)
 description "find the optimal point using golden-search optimization method";
 local R, d, xopt, x1, x2, f1, f2, iter, ea, xint;
 R := (sqrt(5)-1)/2;
 d := R*(xu-xl);
 x1 := xl + d;
 x2 := xu - d;
 f1 := evalf(eval(f, x = x1));
 f2 := evalf(eval(f, x = x2));
 # declare iterator and ea
 iter := 1;
 ea := 1.00;
 # start while
 while ea*100 > es100 and iter < maxiter do
     d := R*d; # new golden ratio
     xint := xu - xl;
     if f1 > f2 then
         xopt := x1;
         fx := f1;
         x1 := x2;
         x2 := x1;
         x1 := x1 + d;
         f2 := f1;
         f1 := evalf(eval(f, x = x1));
     else
         xopt := x2;
         fx := f2;
         xu := x1;
         x1 := x2;
         x2 := xu - d;
         f1 := f2;
         f2 := evalf(eval(f, x = x2));
     end if;
     # calculate new ea
     if xopt <> 0 then
         ea := (1 - R)* abs(xint/xopt);
     end if;
     iter = iter + 1;
 end do;
 # return xopt and fx here
 xopt;
 end proc:

Typing ?while will bring up the help for the Repetition statement, or search for while in the search box under help/Maple help.

As @vv says, there a probably too many solutions, but Groebner gives a way of systematically finding them. You would first want to make all the coefficients rational. As a toy example that finds four real solutions for two equations see below. If you eliminate ones that are unphysical as you go (perhaps negative solutions?), you might be able to find your solution relatively efficiently.

Groebner.mw

(I'm having a problem pasting the contents of the worksheet)

At t=0, you have x(0)=0 and D(x))(0) = 500. Substituting these into eqn gives 0 = 2*500^2, which can't be true. So the equation and the initial conditions are inconsistent at the initial point.

for example, for

f:=(x,y)->x^2;

op(1,eval(f));
gives the sequence of parameters x,y and nops can then be used to count them:

nops([op(1,eval(f))]);

gives 2