@mmcdara You can use this version of BruteForce and then 5/121 works. It also includes your first proc and obtains other representations too, e.g. BF(13/12, 2) = 1/2+1/3+1/4

BF := proc(a, mindenom:=1)
  local A:=a, b:=mindenom-1, s:=NULL:
  while A > 0 do
    b := max(ceil(1/A),b+1):
    A := A-1/b:
    s := s, 1/b:
  end do:
  return s;
end proc:

@student_md 

It's not difficult to see yourself that  int(Psi(s), s=0..t) = P(1)*Psi(t)

implies Psi(t)=0. [it reduces to a simple linear ODE].

@student_md 

Let's take for simplicity the dimension N=1.
Then the first line in Eq (1) ==> Psi(t) = 0. So, your method does not seem to be useful.

@pjv70 

It is hard without sources. If you just want the correct series, simply add 1, or use:

JT3:=convert(JacobiTheta3(z,q), Sum, dummy = n) assuming abs(q)<1:
series(JT3, z=0, 4);

@Carl Love Maple gives only truncated series, not the full ones. It is easy to write the recurrence equations for the coefficients, but rsolve will not be able to solve the system.

@emendes 

It's not clear for me how are you going to use the solutions with the symmetric ones removed. If the system is not symmetric, it will be impossible to know whether (x=a,y=b)  implies (x=b,y=a) too. Also, if the system is symmetric, it is better to treat it as such, returning only (x+y, x*y).

@Kitonum Yes, Mathematica's answer is much better, but my comment was about your remarks which were both incorrect.

@Kitonum 

Not quite (for minimize/maximize):

Edit. The result of Optimization:-Minimize(f, x=0..6) is not a bug, x=2 is a local minimum.
For a global minimum use:
Optimization:-Minimize(f, x=0..6, method=branchandbound);
        [-2., [x = 6.]]

@MapleEnthusiast 

Download puzzle4-vv.mw

@MapleEnthusiast I mean a mathematical reason, not heuristics. Note also that you have simplified wrt size.

(1+e)^3 = 1 + 3e + 3e^2 + e^3: e appears once in LHS and 3 times in RHS. Is this "puzzling"? 

@MapleEnthusiast So, the solutions for A and B are different (and you also took only the numerators!). There is no reason to obtain similar results. If you construct a simpler system (for yourself!) you will be able to check by hand and understand & see the problem.

@MapleEnthusiast Why should they be the same? Note that the R[1,1] could be different in A and B.
If you want to understand exactly what happens, take a smaller system with 2 or 3 parameters and don't use (so many) indices.

@Carl Love Actually, for an expr having a type other than rtable,table, list, Grid:-Map returns 
op(0,expr)(result);
It could be indeed an oversight. Or maybe they want to avoid a (possibly huge) set in the answer.

@acer Yes, I know, but it is anyway unstable; mathematically Rank : R^(n^2) --> N is discontinuous.

So, you are asking for more than the mentioned NP-hard problem!