DEtools:-phaseportrait is restricted to rectangular coordinates. You can get around that by using plottools:-transform to explicitly transform the cordinates to polar after it is passed the plot structure from phaseportrait. This causes some distortion of the arrow shapes, but it's not bad. Finally, the plot structure is wrapped in plots:-display so that we can get polar axes.

If you send the trajectories around the circle again by using t= 0..4*Pi (or more), you'll see that the magenta curve (the outermost) is a limit cycle for all of them---like a lopsided cardioid. But the arrows for the 2nd and subsequent revolutions are not the same, so the plot gets cluttered.

In the code, I use t instead of theta simply to reduce the length of the code.

plots:-display( #needed to get polar axes
    plottools:-transform((t,r)-> r*~[cos,sin](t))( #transform to polar
        DEtools:-phaseportrait(
            [diff(r(t), t) = r(t)*(1-r(t)^2) + 2*r(t)*cos(t)], 
            [r(t)], t= 0..2*Pi, 
            `[]`~(r(0)=~.25*~[$1..4]), #inits
            linecolor= [blue, DarkGreen, black, magenta],
            numpoints= 1000
        )
    ), 
    axiscoordinates= polar, size= [1000$2], thickness= 0
);

Here are two more ways:

fsolve(Statistics:-CDF(Normal(mu,10), 500) = 0.05, mu= 500..infinity);
                          516.4485363

fsolve(Statistics:-Quantile(Normal(mu,10), .05) = 500);
                          516.4485363

In the worksheet below, I show how to get the symbolic solution for part d (c = 2.5) by eigenvalue methods. You'd still need to fill numerous of the step-by-step details.
 

Download MatODEs.mw

There is a Maple command to compute an interpolating polynomial:

''f''~((V:= [6,7,-1])) =~ map[3](CurveFitting:-PolynomialInterpolation, [$1..5], 1/~[$1..5], V);

This uses Maple 2018 syntax.

Regarding your handle, "lovemaple": You should show your love by writing some Maple code.

To avoid having a series of Questions with the same title (which is confusing), I changed you titles.

There is a typo in your Poisson probability formula. It should be P(k) = lambda^k*exp(-lambda)/k!.

Both parts of your problem can be done efficiently with a single 4-line procedure that has a combined for-while-loop (thus combining the exercise's requirements for each part), because both parts require summing the above formula for k = 0, 1, 2, ..., N for either a given N (part a) or for an N to be determined by the sum exceeding a certain value (part b). The efficiency is due to not computing lambda^k or k! explicitly; rather, their quotient is built termwise. 

I urge you to deconstruct the 4 main lines (i.e., you can totally ignore the lines that I commented (in the code) as being for advanced users) of the main procedure below, to figure out what each part does (from help pages or asking here), and to reconstruct your own solution from your new-found knowledge. You may also need a brief tutorial or refresher on discrete probability distributions. Wikipedia is an excellent source, and you can also ask here, even if the question has little-to-none Maple content.
  
This code uses Maple 2019 syntax!
 

Download Poisson.mw

In the worksheet below, your 4th proposed solution violates your 5th equation for any m other than -1, 0, or 1. By direct substitution using m=2, p=-1, q=1, it is explicitly shown to violate the 5th equation.

Your 3rd proposed solution is a specialization (to your proposed value of a1) of the 2nd solution returned by solve, which leaves a1 arbitrary.

Your 1st and 2nd proposed solutions, if you rationalize their denominators, can each be seen to be equivalent to a solution returned by solve.
 

Download RatDenom.mw

It's very difficult to measure the time of such trivial operations because the time of the measuring and looping code overwhelms what's being measured.

My first guess is that nops(data) <> 0 is the fastest way. My close second guess would be

not (data::[]).

Here are some hints:

1. There are four ways that may occur to you:

1. Generate the fibonaccis and test if they're square.
2. Generate the squares and test if they're fibonacci.
3. Generate both and find the intersection.
4. Find a formula that just gives square fibonaccis.

Of these, 1 and 3 are easy, and 1 is the most computationally efficient. 2 is more difficult because it's much easier to check if a number is square than to check if it's fibonacci, and also because there are far more squares than fibonaccis. 4, if it were possible at all, would require some sophisticated mathematical work.

2. The predicate issqr tests whether an integer n is a square, i.e., issqr(n) returns true or false depending on whether n is square.

3. A while-loop can be used to stop generation when you exceed 10000.

4. In Maple 2019, the whole thing can be done in a single line of code.

Use

plots:-spacecurve([0, y, y^2],  y= -1..1);

You're absolutely correct that the series is divergent. In some cases, Maple uses something called "formal summation" to evaluate a divergent series. This is especially likely to happen when you use evalf on an a sum that returns unevaluated or on an inert Sum. You can prevent this by setting

_EnvFormal:= false:

anywhere in your session or in an initialization file. 

The expand command with split the integrals. You could do this:

Int(exp(-I*(A*x+B*y)), [x= -a...a,  y= -a..a]):
convert(value(expand(%)), trig);
 

I think that the HSV color scale provides a more-intuitive way than RGB of translating reals from bounded intervals to colors. The below code is a variation of Acer's, but I scale along H from 0 to .85 and keep S and V constant at 1. H values higher than about .85 look closer to red than violet to me. 

I have also provided what I think is a more intuitive way to use value to get the numeric values.
 

Download TypesetHSV.mw

It can be done like this:

rangetolist:= (x::range(realcons))-> [seq(x)]:

If you want to interpret the instructions literally, replace range(realcons) with (integer..integer).

It might be possible if the expressions are polynomials. See ?solve,parametric. The algorithm used by this command is extremely complicated, and results often take a long time to compute.

Using RealDomain seems to interfere with the ability of simplify to process assumptions. So, don't use it. I don't know exactly why it interferes, but I do know that RealDomain alters the definitions of both sqrt and simplify.

AFAIK, testeq is a probabalistic test that uses modular arithmetic. Because of the modularity, I don't think that it can ever process assumptions. You should replace testeq with is, which is a deterministic test that uses ordinary complex-number arithmetic and is the main command for deriving true/false conclusions from assumptions.