As I indicated to the OP on the other thread, if fsolve doesn't find an answer, it returns NULL (as in the response by @tomleslie) or unevaluated (as in the response by @mmcdara). If you want to test for both these possibilities, then you want

ans:=fsolve(x=x+1,x);

if ans=NULL or op(0,ans)=`fsolve` then "no" end if;

Note this handles also the multivariable case, where a valid solution won't just be a number or sequence of numbers, but something like {x=1, y=3}.

Testing NULL output from fsolve should work here, though sometimes fsolve just returns unevaluated.

Seems to me it is working, just fsolve is really slow. You can speed it up by using the last solution as the initial estimate (but convergence of your outer loop is still really slow).

Try

Q[h+1] := fsolve(eval(f[3], {R = R[h], S = S[h]}) = 0, Q = Q[h]);
if Q[h+1] = NULL then print("breakQ"); break end if;

If you work out Theta[k] and Phi[k] for k=2 nd 3, then you can use a single loop from 4 to M for Theta[k], Phi[k], F[k]. (I shouldn't really cut and paste to duplicate lines here.)

Solution_1-dharr.mw

as @mmcdara pointed out, you van specify ranges in which to search for solutions. You can also use the avoid option to not get the same solution, e.g. after ans1:=fsolve({f, g}); you can use fsolve({f, g},{x,y},avoid=ans1); to get the second solution.

plots:-complexplot(exp(I*theta),theta=0..2*Pi);

put the plot in a variable, e.g., p1:=plot(x^2,x=0..3):

and then lprint(p1);

Even though the first two seem to be interpreted as 2-D integrals, the correct syntax for 2D integrals needs the variables and ranges in a list (or nested 1D integrals), like

simplify(Int(integrand, [phi = 0 .. 2*Pi, r = 0 .. 2]));

and then you get the correct answers. 

You have two first order equations, so if you know omega[0] and sigma, you need only two initial conditions. Put them as a sequence not as a list to solve the syntax error. So give omega[0] and sigma values then

ICS:=C(0)=0,A(0)=1:
sol:=dsolve({C_t,A_t,ICS},{C(t),A(t)},type=numeric):
pRange:=0..20:
plots:-odeplot( sol, [t, A(t)],t=pRange, numpoints=10000 );  #note A(t) not A

If you want to find the values of omega[0] and sigma by specifying two additional conditions, then I think that only works with boundary value problems, and you have only initial conditions. Note also the syntax for derivative conditions

is D(C)(0)=0 for derivative zero at time zero 

I think dchange existed in Maple V, so try

integral := Int(2*(sin(theta)/cos(theta))^(2*p-1), theta = 0 .. (1/2)*Pi);
PDEtools[dchange]({theta=arctan(t)},integral);
value(%);

gives Pi*csc(Pi*p) (in Maple 2017)

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