@torabi Yes, Maple has no numeric solver for elliptic PDEs. You're out of luck.

@John Fredsted I think that it's the best option so far (for a small integer like 7). I didn't hold out any hope of the combinat:-randcomb being efficient.

Is f3(x) a known function? Can h3 be given a value before solving? On a quick read through, I don't see anything that makes this system nonlinear.

@Carl Love Here's a novel approach, and it's a short one-liner.

GenXY:= ()-> combinat:-randcomb(7, 2)[];

This generates a uniform distribution.

@Kitonum The distribution that you're sampling has the same problem as John's. See the histograms that I posted in a Reply to his Answer.

@John Fredsted 

The difference in the distributions being sampled is profound, as the following histograms show. Figuring out why this difference occurs might be a good problem for a junior- or senior-level probability course.

Download two_distributions.mw

@yongling93 In that case, you want to store your bits in a Matrix.  Then all your shifting (what you call "crossover") can be done with surprisingly simple indexing commands. Like this:

B:= Matrix([[0,0,1,1,0,0,1], [0,0,1,1,1,0,0],[0,1,0,1,0,1,0], [1,0,0,1,1,1,0]]):
Bnew:= copy(B):
Bnew[[4,1,2,3], 3..5]:= B[.., 3..5]:
Bnew;

B;

I converted this from a Question to a Post because the question asks for open-ended discussion and opinions rather than definite answers.

@Carl Love Please let me know the outcome of this assignment. I particularly want to know if my guess about the professor trying to stump you with a impossible n was correct.

@adrian5213 Just add ,rn to the return statement.

Do you mean the tensor package? There is no standard package named debever, but ?debever takes you to help about the tensor package.

Here's a Maple worksheet that does it.

Download NH4Cl_density.mw

@adrian5213 Yes, it's easy using the same binary converter that I gave above with a new operator that gives the quotient and remainder.

`&B`:= (ex::uneval)->
     subsindets(eval(subsindets(ex, integer, convert, decimal, 2)), integer, convert, binary):

QR:= proc(N::integer, n::integer) local r, q:= iquo(N,n,r); [q,r] end proc:

&B QR(010100111, 10010);

@tomleslie You don't need to worry about options c[3], c[4], c[5]. The first, c[3], only applies if itask = 4 or 5; you're using itask = 1 (the default). The other two, c[4] and c[5], only apply if you're using one of the -band submethods of lsode; you're using submethod adamsfull. These options specify how many super- and sub-diagonals, respectively, of the jacobian to use. Personally, I think that using a jacobian-based method for such a "wild" RHS is a bad idea. I think that the submethod should be changed to adamsfunc, which uses no jacobian.

@Kitonum Sigh. I guess that it's yet another difference between 1D and 2D input. I guess that as the 2D input is internally translated to 1D, an extra level of evaluation occurs.