To get all 50 roots:

degree(g,x); #50
fsolve(g=0,x,complex); #Include the imaginary ones
nops([%]); #50

@canoiras Although it is certainly not clear to me precisely what you want to do, I should point out that in the example, which I gave in my answer below, I used boolean variables essentially to switch from one form of the system to another. These boolean variables are changed by the events specified.
Maybe you could give the Maple code, you already have, for us to work on.

@Bendesarts You cannot have equations in plot. Use a parametric plot:

plot([seq([Pz[i],t,t=-0.1..0.1],i=1..nops(Pz))],colour=yellow,linestyle=3, thickness=2,view=[0..2,default]);

Adjust [Pz[i],t,t=-0.1..0.1] to taste.

@Bendesarts If you do this:

restart;
for i to 4 do
  u[i]:=unapply(L/2*cos(w*t+phi[i]),t);
  Equ[i]:=evaln(u[i](t))=u[i](t)
end do;
## then you have
Equ[1];
              (1/2)*L*cos(t*w+phi[1]) = (1/2)*L*cos(t*w+phi[1])

Certainly true, but not exciting!

So if you just want to define the u[i]'s as functions, then you shouldn't define equations, you should just assign to u[i] as in e.g:
restart;
for i to 4 do u[i]:=unapply(L/2*cos(w*t+phi[i]),t) end do;
#Test:
u[1](s);

@Thomas Richard This way it exhibits well as it is run, but try afterwards:
Eq[1];

This is just stating that a=a, so to speak.

I tried from within Maple using "Check for Updates". It took forever (50 minutes at least), and stalled when 8 seconds were left, or so it claimed. Then at some point it "decided" to try again and announced a long download time (about 50 minutes again). I went shopping for beer. When I came back my computer had fallen asleep, but when I woke it up the downloading was not finished, but later it stalled at 5 seconds left. Then apparently it gave up. 

I then tried the Maplesoft website. The update downloaded very rapidly and installed without problem.

There appear two versions of _Z1 and two versions of _Z2, and that is weird. Certainly a bug.

restart;
eq:=Diff(W(t), t) = -q*V*(sin(Phi)-sin(Psi[s]))/(2*h*Pi)
soln := solve([rhs(eq) = 0, Phi < 2*Pi], [Phi], allsolutions = true,explicit);
idts:=indets(soln,`local`);
nops(idts); #4
addressof~(idts);
nops(%); #4
getassumptions(idts[1]);#
getassumptions(idts[2]);#These are not the same assumptions
Lidts:=convert(idts,list);
S:=Lidts=~[z11,z12,2,3];
subs(S,soln);

Note: Another strange thing: If you replace Psi[s] with an unindexed name like psi, then you will notice that idts has five members: assumptions are made on psi. If you try replacing Psi[s] with Psi (which is an initially known function in Maple) you get an error message from solve saying
Error, (in assume) cannot assume on an assigned name

Please see a related discussion:
http://mapleprimes.com/questions/208328-Differential-Equation-And--Notation-Question

Get rid of the units (after having made them conform to each other) or introduce dimensionless variables. Then compute.

@tomleslie I'm not a big fan of units either. 
The expression exp(-4000*t) can be nothing but a pure number with no unit no matter what t is or what unit is used.
Just think of the series expansion for exp.
exp(x) = convert(exp(x),FormalPowerSeries,x);

If the unit of x is U what is the unit of exp(x)?

My thinking is: get your units straight from the start, perhaps with the help of Maple, then forget about them, and compute!

@tomleslie What would be the unit for t?

@shahid Sorry, I don't.

@shahid Since you have a Maple version earlier than Maple 18 eval[recurse] doesn't work.

This should work (I did this in Maple 15):

restart;
eq1:=u(x,y)*(diff(u(x,y), x))+v(x,y)*(diff(u(x,y), y)) = (diff(u(x,y), y, y))/(1+epsilon*theta(x,y))-epsilon*(diff(u(x,y), y))*(diff(theta(x,y), y))/(1+epsilon*theta(x,y))^2+theta(x,y);

eq2:=u(x,y)*(diff(theta(x,y), x))+v(x,y)*(diff(theta(x,y), y)) = (diff(theta(x,y), y, y))/Pr;

## Replace the original S with the following:
S1:={u(x,y) = diff(psi(x,y), y), v(x,y) = -(diff(psi(x,y), x))};
S2:={ psi(x,y) = x^(4/5)*f(x,eta)/(1+x)^(1/20), theta(x,y)=x^(1/5)*(1+x)^(-1/5)*h(x,y)};
S3:= eta = y/(x^(1/5)*(1+x)^(1/20));
S23:=eval(S2,S3);
eval([eq1,eq2],S1);
eval(%,S23);
subs(rhs(S3)=lhs(S3),%);
collect((lhs-rhs)~(%)=~0,[D,diff],factor);
res:=subs(y/(x^(1/5)*(1+x)^(1/20))=eta,%);

##To make it more readable you can continue with:
PDEtools:-declare(f(x,eta),h(x,y));
convert(res,diff);

About rewriting S1*S2-S3*S4.

You can rewrite the products of the inert sums S1*S2 and S3*S4 using the Cauchy product. Thus you get

S1*S2-S3*S4 = S12 - S34, where S12 is the Cauchy product of S1 and S2, and S34 is the Cauchy product of S3 and S4.

The difference S12-S34 can be combined to just one sum S:

S := 1/(2*Pi)*Sum(x^(3*k1+2)*3^(1/2+k1)*(Sum((3*GAMMA(k+2/3)*(-2/3-k1+k)*GAMMA(k1+1/3-k)+3*GAMMA(4/3+k)*GAMMA(2/3+k1-k))/(GAMMA(3*k+2)*GAMMA(3*k1-3*k+3)), k = 0 .. k1)), k1 = 0 .. infinity);

This can be written in short form as
1/(2*Pi)*Sum(x^(3*k1+2)*3^(1/2+k1)*m(k1), k1 = 0 .. infinity);

where m(k1) < 0 for all k1 >=0. It is given by:

m(k1) = Sum((3*GAMMA(k+2/3)*(-2/3-k1+k)*GAMMA(k1+1/3-k)+3*GAMMA(4/3+k)*GAMMA(2/3+k1-k))/(GAMMA(3*k+2)*GAMMA(3*k1-3*k+3)), k = 0 .. k1);

m(k1) is very rapidly decreasing in absolute value.

Now you don't have cancellation problems for x > 0 since all terms in S are negative.
But I wonder why x-values as high as x=100 are interesting to you. S will be huge (but negative) for x = 100.
For x=10 you get
evalf[50](eval(S,x=10));
    -4.7714630011989009659113681348806733279336141633480*10^8
and for x = 11:
evalf[50](eval(S,x=11));
    -1.1891747657696851151848674770800968156039786693327*10^10
### Notice that you will run into similar problems with the well known power series for exp(x):
restart;
S:=Sum(x^n/n!,n=0..infinity);
seq(evalf[500](eval(S,x=i*10)),i=0..10);

@shahid Here is what I got (last output shown only) when I just now copied my own code from above and executed the whole thing in Maple 2015.2. Actually, in my worksheet it is shown with the abbreviated notation used because of my use of PDEtools:-declare(f(x,eta),h(x,y)). But when I copy and paste it into MaplePrimes it shows up as you see below.

You could begin by telling us what Ni1 is. Maybe you also can tell us if a1 and a2 are assigned to anything.