@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.

@Markiyan Hirnyk Maybe it's a PDF if x is restricted to x >= 0. Many (perhaps most) distributions have such a restriction.

@JohnS The scanf family of commands (fscanf, sscanf, and scanf) always return a list. Even if all the scans fail, they return the empty list; if you're scanning only one item, they return a list with that one item.

@Josolumoh Hold your horses. You've been very slow at providing the necessary details. Is it a continuous or discrete distribution? What is the support, i.e., the domain for x? What are the domains of the parameters b, d, and r?

The appropriate command is Statistics:-Sample.

@bloodtalon What is your Maple version? Like I said, it won't work with Maple earlier than 18, when the key option to sort was introduced. But Markiyan's adjustment (immediately above)  will make it work for earlier Maple.

@Josolumoh Is your function a probability density function (pdf)? If so, what's the name? If it's a standard one, Maple likely already has a sampling procedure for it.

@JohnS If L is a list or set, then L[] extracts the underlying sequence. This is documented in the ninth paragraph of Description at ?list, which simply states

To extract the contents of a list or set, use the empty selection operator [ ].

That's not possible: The function is much too complicated. Perhaps if you describe why you want to invert it, we could come up with an alternative.