@mmcdara V=~0 is not a solution because of 0/0, but perhaps it is a solution to the original problem before it was manipulated into this form. I agree there is something strange about the system; perhaps for other values of the OP's free parameters it has clearer solutions.

@mmcdara I dont understand your final comment, since your minima have mu__jk=mu__ki; they also have the lambdas with the same magnitude. You are right this is a tricky problem for fsolve (so needs careful checking increasing Digits and perhaps the fulldigits option), but minimizing only finds one of the solutions, and prseumably all solutions are sought.

If you put all the equations in numer/denom form with normal, then take the numerators, then all variables=zero is a solution:

V=~0;
eval(map(numer@normal,Eqp),%);

gives zeroes. The numerators are all multivariate polynomials, though not sure if this helps.

@MaPal93 Here's how I would collect output of fsolve into a Matrix: fsolvecollect.mw

I think your premise that the different solutions of fsolve differ because of convergence issues is not correct; they are separate accurate solutions of the set of equations. To verify this, you could set Digits:=20; and see that almost all of the digits you currently get are accurate.

To put values back into alpha etc use eval: eval(alpha,{mu=2,s=3})

To vary the parameters means putting things into some sort of loop; a lot depends on how systematic you want to be.

@MaPal93 It seems very unlikely there is an analytical solution; solve already will try many algorithms.

@MaPal93 

• what the ^+ does at the SMI step (4th line)
  It means transpose (%T is an alternative), so SMI is a row vector, as you wanted.
• why in natlq2 and natlq3 steps SMI2 and SMI3 are not transposed (see your initial comment) before being multiplied by Y
  I think the above probably already answers that. The row vectors SMI2 and SMI3 multiply the column vector Y to give a scalar. (The usual dive row down column and add, which in this case amounts to a dot product.)

For fsolve, you get one root (at least in my earlier version of Maple). To find others you will need to help fsolve. If you are lucky, use the avoid option to avoid this one next time and then so on. Or you can give ranges within which you can find the root of interest, or you can give initial guesses (less reliable in my experience). Maple has optimizers if you need to go that route.)

See my answer here. But also maxsols=8 may work without these tricks.

@MaPal93 I don't understand; your newly updated document has no LinearSolve. I used worksheet mode originally because I am most familiar with it, but have done the updated one in document mode - see the added document to my answer above. Run the document, and remove the colons to see the answers. The answers are extremely long; I'm not sure they will be useful.

ans1 is a 1..3, 1..13 Array, as may be seen by changing the : to ;. So ans1[12..14,3] leads to an indexing error, since 12 is not in the range 1..3.

By the way, why not just export the Array contents directly, rather than extracting them from a plot?

@ijuptilk Sorry, forgot about discont only for 2D. But at least in Maple 2015 I am not seeing a problem, though I didn't explore all ranges.

plot3d(y/(x+1)+1/(y^2-y), x = -3 .. 3, y = -3 .. 3, view = [default, default, -8 .. 8],
 style = surfacecontour);

gives

@ijuptilk use discont=true if there are discontinuities.

@Carl Love Excellent effort, going beyond. Vote up.

@mehdibgh Here's an example - words refers to memory usage. (The second PrintProfiles is for both calls to test.)
 

Download Profile.mw

Maple has the CodeTools:-Profiling package to help you find where these limitations are. As a first step, you can put the critical part in a procedure and pass it to CodeTools:-Profiling:-Profile - see the help page for that.

@Ronan For URLs I think you need forward slashes always. For Maple filenames, forward slashes always work (platform independently), but on windows backslashes work, but need to be escaped, i.e., \\ is \

(I said the text file is opened by notepad, but actually it is by the browser.)

@jeremy_murphy Maple's combstruct package implements much of the material in Flajolet's book. One defines a combinatorial system, e.g., binary trees, using a (usually recursive) grammar, and then you can calculate the number of objects (trees) of a given size n (=number of leafs), list all objects of a given size, and find a generating function whose coefficients are the number of objects of the different sizes. The entropy would normally be ln(number of objects) (this is the Shannon formula when all objects have the same probability), or sometimes one has a "entropy constant" (see "hard hexagon entropy constant" on mathworld,com, which is the limit((number of objects of size n)^(1/n), n->infinity), or exp(entropy/size).

An example such as sequences of 1's and 2's that sum to n might be closer to what you are thinking, since there could be different probablities for sequences with different numbers of 2's.

The combstruct package is quite powerful but needs some practice to figure out. If you can be more specific about what you want, or want me to show an example I can do that. I think, however, that getting from the generating function to the entropy is non-trivial.

@Carl Love Number of states is clearer. I was confused about the use of n for the size of the combinatorial "object", and for the size of the "outcome alphabet". In the sorts of physics objects I tend to deal with, the size could be say n=5 for 5 sites, and then the number of states is the number of possibilities (with equal probabibilities or not) and is larger than that. The formula for entropy is the same. For equal probabilities the formula simplifies to ln(W), W = number of possibilities. The derived formula I gave is then the limiting entropy per site for n infinity. In this case the result can be viewed differently as a two-state system as you pointed out (we have moved from a micocanonical ensemble picture to a grand canonical ensemble picture in physics terminology).

Perhaps the OP can give some more specifics about the type of functions. But I'm guessing if any software can find a symbolic solution, Maple can!