The number of which primes?  Do you mean something like this?

for a to 10 do
  for b to 10 do
     for c to 10 do
        printf("There are %d primes of the form %a for x from 1 to 100\n", 
           nops(select(isprime, [seq](a*x^2+b*x+c, x=1..100))), a*x^2+b*x+c)
end do end do end do:

The problem is that the result of inequal consists of a background in the "feasible" colour, overlaid wit  h polygons in the "excluded" colour for each inequality.  The "feasible" colour shows through in the region that's not covered by any of the "excluded" polygons.  But combining "inequal" plots with display will only ever show you one of the plots; the others will be hidden.

If you want to see feasible regions from several sets of inequalities, you'll need a different approach that will colour in only the feasible region.  This can be done using implicitplot.  For example:

> with(plots):
  implicitplot(max(x^2+y^2-1, x - y) <= 0, x=-1..1,y=-1..1,
      coloring=[red,white],filledregions, gridrefine=3,crossingrefine=3,
      transparency=1/2); P1:= %:
  implicitplot(max(1/2-x^2-y^2, x*y - 1/4) <= 0, x=-1..1,y=-1..1,
      coloring=[blue,white],filledregions, gridrefine=3,crossingrefine=3,
      transparency=1/2); P2:= %:
  display([P1,P2]);

1) I don't know what you mean by "floor doesn't take into account decimals".  It does.

2) Hint: fib(n) is very close to phi^n/sqrt(5) where phi = (1+sqrt(5))/2.

Note that x^2 >= x^3 on the interval in question.

> plottools[transform]((x,y) -> [x,y+x^3])(plot(x^2-x^3, x=0..1, filled=true));

This will work in the Standard GUI, but not Classic. 

Or, for one that works in Classic as well:

> plots[implicitplot](max(y-x^2,x^3-y), x=0..1, y=0..1, filledregions=true,
    coloring=[red,white], grid=[50,50], gridrefine=3, crossingrefine=2);

You got the series B, but you didn't use it in C, you used D(f)(x) instead.  I think you want

> N:= series(x - A/B, x=a, 3);

N := a+1/2*`@@`(D,2)(f)(a)/D(f)(a)*(x-a)^2+O((x-a)^3)

And then

> e[n+1] = subs(x-a = e[n], series(N - a,x=a,3));

e[n+1] = 1/2*`@@`(D,2)(f)(a)/D(f)(a)*e[n]^2+O(e[n]^3)

First of all, you need to be clear on what region you're talking about.  There is no "region bounded by y=x^2 and line x=1".  I suspect what you mean is the region bounded by y=x^2, y=0 and x = 1.  For that you could do

> plot(x^2, x=0..1, filled=true);

"Damped, overdamped etc." are properties of a linear system.  In general with a nonlinear system, you may have several critical points (equilibrium solutions), which may differ in their properties.  In your case, since the exponential function is never 0, there are no critical points at all.  I don't know what sort of classification you're looking for in a case such as this.
 

For example:

with(Maplets[Elements]):
maplet := Maplet(FileDialog['FD1'](
    'filefilter' = "txt",
    'filterdescription' = "Text Files",
   'onapprove' = Shutdown(['FD1']),
   'oncancel' = Shutdown()
)):
Maplets[Display](maplet);

rifsimp in the DEtools package will do it, but the result is pretty complicated.  Note: I'm leaving out the x because you have no differentiations with respect to x.

> sys:= [P(y)+P(y)^3 + Q(y)^2 - R(y) = 0, Q(y) + Q(y)^2 + P(y)^3 + S(y)=0,
         diff(P(y),y)+S(y)=0,diff(Q(y),y)-R(y)=0];

> DEtools[rifsimp](sys,[P(y),Q(y)]);

Well, I see you've changed some coefficients...  This time to avoid the singularity in the polylog you need y < 23.713 approximately, so at least it's not so obviously inconsistent with the boundary condition at infinity.  But it still won't work.  If y(r) -> 0 as r -> infinity, the DE becomes approximately

diff(y(r), r, r) = -2*(diff(y(r), r))/r −363801.6496

and all solutions of this have y(r) -> -infinity as r -> infinity.

It would be full evaluation, but the quotes in 'a' prevent the evaluation from going further.

In a case such as this, where the results are naturally organized as a two-dimensional array of points (one dimension for s, the other for tau), I think it's better to use surfdata in the plots package, rather than SurfacePlot which needs to figure out by itself how to organize the data.   There are also several other issues:
- it's extremely inefficient to produce long lists the way you're doing it, appending one element at a time.
- use add rather than sum for adding a finite number of objects
- do as little as possible inside an inner loop: in particular avoid defining functions there
- use hardware floats whenever possible. 

You could try something like this:

> A:= Array(1..50, 1 .. 7201, 1..3, datatype=float[8]);
   for i from 1 to 7201 do
      tau:= (i-1)*0.025;
      for s from 1 to 50 do
         A[s,i,1]:= s;
         A[s,i,2]:= tau;
         A[s,i,3]:= evalhf(add(0.25e-2*f[m+1]*(1-(0.25e-2*m-tau)^2/s^2)*exp(-(1/2)*(0.25e-2*m-tau)^2/s^2)
             /(sqrt(2*Pi)*s^2), m = 0 .. 71999));
    end do end do:
    with(plots):
    surfdata(A, style=patchnogrid, axes=box);

But there's still a problem here: you'll have 50*7201*72000 = 25923600000 individual terms to add up. 
Even on a very fast computer, that will take some time. Perhaps you could be a bit less ambitious, to reduce these numbers
somewhat.

You might also look at the DiscreteTransforms package, which does wavelet transforms.

You might do something like this.  Given two points A and B, the following procedure breaks up a list L of points wherever a given segment A,B occurs.

> breaker:= proc(L, A, B)
        local R, n,m, breakpoints;
        R:= NULL;
        n:= nops(L);
        breakpoints:= select(i -> (L[i]=A and L[i+1]=B), [$1..n-1]);
        m:= nops(breakpoints);
        if m = 0 then L
        else L[1..breakpoints[1]], seq(L[breakpoints[i-1]+1..breakpoints[i]], i=2..m),
              L[breakpoints[m]+1..n]
        fi
    end proc;
                

Now to remove any segments from point A to point B from a plot structure P, this should work:

> evalindets(P, listlist, breaker, A, B), 

Note that A and B should be exactly in the form that appears in the plot structure.  For example:

> P:= plots[pointplot3d]([seq([i,cos(i*Pi/10),sin(i*Pi/10)],i=0..40)],colour=blue,style=line):
  evalindets(P, listlist, breaker, evalf([21,cos(21*Pi/10),sin(21*Pi/10)]), 
                                   evalf([22,cos(22*Pi/10),sin(22*Pi/10)]));

I assume you're using the version of the geometric distribution where 0 is not allowed, i.e. the discrete distribution

Probability(X = x) = p*(1-p)^(x-1) for integers x >= 1. 

The fact that it's a discrete distribution complicates matters.  For the ratio Z = X/Y, where X and Y are independent random variables with this distribution, the possible values are all positive rational numbers.  If m and n are coprime positive integers, Z = m/n whenever X = m*k and Y=n*k for some positive integer k.  So here is Probability(Z = m/n):

> Px:= x -> p*(1-p)^(x-1);
   sum(Px(k*m)*Px(k*n),k=1..infinity) assuming p>0,p<1,m::posint,n::posint;

p^2*(1-p)^(m-2+n)/(1-(1-p)^(m+n))

The CDF will not have a closed form, but can be done numerically: it is

> F:= x -> evalf(Sum(floor(x*s/(1+x))*p^2*(1-p)^(s-2)/(1-(1-p)^s), s=2..infinity));

For reasonable values of p, we can get a very good approximation by using a partial sum, say for s = 2..1000.

> Fapprox:= subs(infinity=1000, eval(F));

> p:= 1/2; plot(Fapprox, 0 .. Pi);

The reason you don't have arrows is that your system is not autonomous: it depends explicitly on t.  The arrows are supposed to show you the direction trajectories will go, but this changes with time.