@AfshinK I figured it out anyway. Thanks.

Lest anyone think that the summation or the subscripts are the source of your error, consider the simpler case c=1, and remove the subscripts:

solve(m*exp(-(m-1)*t)*d^2-theta+theta*t-theta*t-m*eta*(1-exp(-t))^(m-1) = 0, t);

This produces the same error message, so the summation and the subscripts have nothing to do with it.

@Christian Wolinski By "remainder" I assume that you mean residual. I wouldn't expect to get any better residuals from factor than from fsolve. But the fsolve residuals aren't bad (in Maple 16):

Digits:= 500:
Sols500:= [fsolve(rho_poly, complex)]:
unapply(abs(rho_poly), rho)~(Sols500):
max(evalf[5]~(%));

Also consider the absolute errors:

Digits:= 600:
Sols600:= [fsolve(rho_poly, complex)]:
abs~(Sols600 -~ Sols500):
evalf[5](max(%));

And the relative errors:

abs~((Sols600 -~ Sols500) /~ Sols600):
evalf[5](max(%));

What does the above code return in your Maple?

@acer Your first two examples contradict each other: The input is the same, but the output differs. I think that you may have made a copy-and-paste error.

Why use factor? You can use fsolve(rho_poly, complex) to get all the roots.

@Markiyan Hirnyk I don't understand why one way works with IsPrime and another way with PrimeDecomposition. I suggest going with Axel's way.

@bloodtalon Directly above, you say that you only need the x*y*z values. Earlier, Markiyan gave you a command that gives just the extremal values: extrema. To that you said (I'm quoting directly)

I need to know the points where these are achieved

So, you contradict yourself.

Maple 2015 is the latest version. Maple 18 was the version before that.

@Markiyan Hirnyk 

IsPrime(J, sqrt2);

and

IsPrime(J, sqrt(2));

both return false.

Have you read the help ?PolynomialIdeals,IsPrime? An example very similar to this is shown.

Or, if you want to explicitly declare it in the ideal:

J:= <x^2*y^2 - 2, sqrt2>:
IsPrime(J);

     false

@Josolumoh You could put a coefficient on the whole thing that would make the sum 1. The coefficient would need to computed for each (d,b,r) triple.

@vv Okay, I see. My program was just finding some arbitrary small value of d such that the sum is 1 to ten decimal places, and actually any smaller nonnegative value would do the same. Yes, I see now that only d=0 makes the sum exactly 1 for any b and r.

@nMaple Your inequality Y <= 20 contradicts your code fragment

for (double Y = 0; Y<=10; Y+=0.5){

Which is correct, and what is the step for Y?

I guess that you want to plot six separate surfaces, right? Since you have six separate columns of dependent-variable data, that's the only way that I can make sense of your request.

@Josolumoh I already anticipated this problem, and I had put my Reply to it below the Answer below. Look there.

@Josolumoh Okay, for your function note that the range given in the Sample command needs to be finite even if the Support given in the Distribution command is infinite. In this case, you choose the range such that the probability of being outside the range is insignificantly small.

@patient Your so-called "maximum" is a singularity. You can make it as close to infinity as you want by getting closer to the origin (alpha=0, x=0). But, depending on its error-control settings, dsolve gives an error if you get too close. That is what happened to the plot on the left.