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

@mskalsi Your errors have nothing to do with my code. Your first error, Error, recursive assignment, happened because you tried to use A in two ways. You wanted it to be both the name of the Matrix and one of the coefficients and hence one of the entries of the Matrix. A Matrix can't contain itself as own of its own entries. An attempt to do it is called a recursive assignment.

Your second error is caused by essentially the same thing, but in a more subtle and insidious way. You are trying to use both A[1] and A as separate variables. But they aren't separate, as the following shows:

A[1]:= 7:
A:= 3:
A[1];

     

If you want to use both A and a subscripted form of A in the same program, then use A and A__1 (the 1 can be replaced by anything). The A__1 will appear subscripted in the output (and also in the input since you're using 2D input). These will be treated as completely separate variables.

Regarding the disordering caused by indets: The order in the final Matrix is determined by the order in the Vector product <U[1],U[2]>.<a[1] | a[2]>; the order returned by indets is irrelevant.