A correlation coefficient does not determine the joint distribution of the coins.  I'll assume that this distribution is symmetric with respect to permutations of the coins.  Then there are three degrees of freedom: we may take P(HHH) = a, P(HHT) = P(HTH) = P(THH) = b, P(HTT) = P(THT) = P(TTH) = c, P(TTT) = d where a,b,c,d >= 0,
a + 3*b + 3*c + d = 1.  If X_i are the Bernoulli random variables corresponding to heads for coins 1,2,3, then the correlation coefficient of any two of these is
 

corr := (a+b-(a+2*b+c)^2)/(a+2*b+c-(a+2*b+c)^2)

You could specify the correlation coefficient and the expected value (i.e. the probability of heads for each coin), and still have one degree of freedom.
Thus with correlation coefficient p and expected value m:
 

> eqs:= {a+3*b + 3*c + d = 1, a + 2*b + c = m, corr = p};
  solve(eqs,{a,b,c});

{a = -3*m-d+1+3*m^2+3*p*m-3*p*m^2, b = 3*m+d-1-2*m^2-2*p*m+2*p*m^2, c = -2*m-d+1+m^2+p*m-p*m^2}

The probability of three heads or three tails is then a + d = -3*m+1+3*m^2+3*p*m-3*p*m^2.

It doesn't work in Maple 13.02, but here's what I get in Maple 14:

RealRange(Open(-7.203549868),Open(-5.362820746)), RealRange(Open(-4.061957215),Open(-2.221228092)), 
RealRange(Open(-.9203645613),Open(.9203645613)), 
RealRange(Open(2.221228092),Open(4.061957215)), 
RealRange(Open(5.362820746),Open(7.203549868)), 
RealRange(Open(8.504413399),Open(10.34514252))

I don't know why Maple chooses these particular intervals (of course, H is periodic with period Pi, so there are infinitely many intervals) .
 

Basically, I think Maple is not as good at solving inequalities as it is in solving equations.  You might compare:

> solve(H = .9, AllSolutions);

-.9203645613+4.982321776*_B1+6.283185307*_Z1, .9203645613+1.300863531*_B2+6.283185307*_Z2

Of course the solution set for the inequality are some of the intervals with these endpoints.

I'm not quite sure what your Maple code is doing, but I think you may want to do something like this.  Suppose your lists are A and B, each item of each being a list of two items, and you want the list of all items [a[1], a[2], b[1]] where a and b are elements of A and B respectively such that b[2] = a[2].
Then you can use

> map(op,[seq(map(b -> [a[1],a[2],b[1]],select(b -> b[2] = a[2], B)),a=A)]);

Your function decays to 0 extremely rapidly as x increases, because the exponential factor is approximately exp(-559*x).  If you plot this
for x from 0 to 10, nearly all of the surface will be effectively 0.  Try

plot3d(u(x,t), x = 0 .. 0.01, t = 0 .. 10, axes=box, grid=[50,200], style=patchnogrid);

But I suspect something is wrong with your constants.

First of all, you shouldn't use O (which has a special meaning in terms of series, and causes a "too many levels of recursion" error).  Try o instead.

Now it looks to me like there is a real singularity there: f(z)^4*(e(z)^2-.27*(1+z)^3)^2 hits 0.0081 at approximately z = .76663139, after which you'll be taking the square root of a negative quantity.  If you want any chance of a sensible continuation past that singularity, perhaps you'll want to insert an absolute value.  Thus

> eqn10 := diff(e(z)^2, z) = ((3/2)*e(z)^2-3*(n^2/f(z)^2)^2/(2*(e(z)^2-o*(1+z)^(3*g)))
     +(1/2)*((g-1)*3)*o*(1+z)^(3*g))/(2*(1+z));
  eqn11 := diff(f(z), z) = sqrt(abs(1-(n^2/f(z)^2)^2/(e(z)^2-o*(1+z)^(3*g))^2))/(e(z)*(1+z));
  eqn12 := subs(o = .27, g = 1, n = .3, eqn10);
  eqn13 := subs(o = .27, g = 1, n = .3, eqn11);
  sol4 := dsolve([eqn12, eqn13, f(0) = 1, e(0) = 1], numeric,method=classical[rk4],stepsize=0.00001,maxfun=10^6);
  plots[odeplot](sol4,[[z,e(z),colour=red],[z,f(z),colour=blue]],z=0 .. 2,numpoints=1000);

As you can see, something strange happens to f(z). 
It turns out that (e(z)^2-.27*(1+z)^3) hits 0, resulting in df/dz going to infinity, at around z = .7867866833.  Because of the square root, I think the singularity should still be integrable, so f won't go to infinity, but it's going to be very hard for numerical methods to accurately reflect what does happen. 

If you want both spheres at the same time, 

> with(plots): with(plottools):
  display([seq](display([sphere([1,1,1], r), sphere([3,3,3], r)]), r = 0 .. 1.73, 0.173), insequence=true); 

You could fit using trig polynomials.

a := [[1, 1], [3, 2], [3.5, 4], [4, 6], [5, 5], [6, 7], [7, 6], [8, 5], 
[9, 5.5], [10, 4], [11, 1], [12, -5], [11.5, -6], [12, -12], [10, -10], 
[8, -14], [7, -10], [3, -10], [2, -5], [1, -8], [0, 0], [1, 1]];
with(Statistics):
n := nops(a);
Xfit, Yfit:= seq(
  Fit(add(A[i]*cos(2*Pi*i*t/(n-1))+B[i]*sin(2*Pi*i*t/(n-1)),i=0..5),<(j$j=1..n)>,map2(op,m,a),t),
    m=1..2);
plots[display](plot(a,colour=blue),plot([Xfit,Yfit,t=0..n-1]));

From the error message, it looks to me like somewhere in the part of your code that you haven't shown us, you have assigned the value 1 to the variable x.  This doesn't happen (at least in Maple 14) with the code you gave us.  Of course, you need with(Optimization).  You should also use Int rather than int, so Maple doesn't waste time trying to integrate a complicated expression symbolically.  With that, the error message is "invalid arguments".  I think Maple may be confused by the dummy variable of integration: I tried the operator form

> NLPSolve(E,initialpoint=[.1035805, 1.001279, -.561381]);

but got

Warning, no iterations performed as initial point satisfies first-order conditions


                                  [0.103580500000000006 ]
                                  [                     ]
            [7.71495067895969999, [ 1.00127900000000003 ]]
                                  [                     ]
                                  [-0.561381000000000019]

which is rather strange because the gradient vector there is certainly not 0.  I think there may be a bug here.

I'm confused.  Your 1-BIN(z-2,n,p(x,t)) doesn't appear to depend on the integration variable X, so the integral is easy.  Or is X the same as x?
And what equation are you trying to solve, for what variable?

The package name is Statistics, not statistics. 
But your error is premature evaluation: X||i(1,j) won't work with j as a symbolic variable.  The cure is to use add instead of sum.

For example:

with(MultiSeries):
X1 := a1*exp(r1*t);  # growth rate r1
X2:= a2*exp(r2*t);  #  growth rate r2
Avg:= (X1+X2)/2;   # the average
asympt(diff(Avg,t)/Avg, t) assuming r1 > r2;  # asymptotic growth rate of the average

r1+O(1)*(1/exp(t))^(r1-r2)

Yes, the Array arr only contains the z values.  the GRID structure containing that Array has its first and second operands
a..b, c..d, corresponding to the intervals for the x and y variables.  The grid points are equidistant, so if arr is 1..m by
1..n, the data points are [a + (i-1)*(b-a)/(m-1), c + (j-1)*(d-c)/(n-1), arr[i,j]] for i from 1 to m and j from 1 to n.

If you make your plot3d parametric, e.g.

plot3d([x, y, f(x,y)], x = a .. b, y = c .. d)

instead of a GRID structure you will get a MESH containing an Array(1..m, 1..n, 1..3) where the [i,j,1..3] entries are the x,y and z
values for the data points.  You may find this more convenient to work with.

Strange: those equations should be easily solvable.  Perhaps this is a bug.  You can do this:

> D1:= solve(e1,D);
  Qs:= solve(eval(e2,D=D1));
  Ds:= eval(D1,Q=Qs);

The problem may be a form of premature evaluation.  You can try solving an equation involving a procedure rather than an expression.

> q:= unapply(CDF(X,100) - 0.9, n);
  fsolve(q, 0 .. 1);

It seems to be taking a long time...

The right side of the equation X1 is rhs(X1).  However, solve will not work on this equation, which involves both exponentials and trig functions: there is probably no closed-form solution.  What you can do is

> fsolve(rhs(X1)=1, t = 0 .. 1);

(and note by plotting that there is no smaller solution)