This problem is rooted in the combinatorial problem of partitioning an integer n (here N) into at most k parts (here k=3). For details, see the attached entry of Wikipedia as well as Chapter 6 of Ref. [1]. Once you have the list of the partitions of N into at most 3 parts (for N a specified integer, you could use IntegerPartions[N,3] in Mathematica©), using some straightforward combinatorics, you can calculate the expression for the desired energy. Using this expression, you can subsequently solve your specific problem.

For illustration, let us consider the specific case of N=5, and k=3. For the above-mentioned partitions one has [2]:

{5,0,0}, {4,1,0}, {3,2,0}, {3,1,1}, {2,2,1}.

-- {5,0,0} implies that the available 5 particles can be in one of the 3 energy levels. For this there are 3 distinct possibilities: [5,0,0], [0,5,0], [0,0,5]. Here [i,j,k] denotes a state with i particles in level 1, j particles in level 2, and k particles in level 3.

-- {4,1,0} implies that 4 particles can be placed in one energy level and the remaining 1 particle in another energy level. For this there are 6 distinct possibilities: [4,1,0], [4,0,1], [0,4,1], [1,4,0], [1,0,4], [0,1,4].

-- {3,1,1} implies that 3 particles can be placed in one energy level and the remaining 2 particles in two remaining energy levels. For this there are 3 distinct possibilities: [3,1,1], [1,3,1], [1,1,3].

-- {2,2,1} implies that two pairs of particles can be placed in two distinct levels and the remaining one particle in the remaining level. For this there are 3 distinct possibilities: [2,2,1], [2,1,2], [1,2,2].

Denoting the total energy by E, and those of the three levels by E1, E2 and E3, with Ai = exp(-Ei/kBT), i=1,2,3, for the problem at hand one has (under the assumption that particles [here bosons] are indistinguishable):

E = [(5 E1 A1^5 + 5 E2 A2^5 + 5 E3 A3^5)

+(4 E1 + E2) A1^4 A2 + (4 E1 + E3 ) A1^4 A3+ (4 E2 + E3) A2^4 A3

+ (E1 + 4 E2) A1 A2^4 + (E1 + 4 E3) A1 A3^4 + (E2 + 4 E3) A2 A3^4

+(3 E1 + E2 + E3) A1^3 A2 A3 + (E1 + 3 E2 + E3) A1 A2^3 A3

+ (E1 + E2 + 3 E3) A1 A2 A3^3

+(2 E1 + 2 E2 + E3) A1^2 A2^2 A3 + (2 E1 + E2 + 2 E3) A1^2 A2 A3^2

+ (E1 + 2 E2 + 2 E3) A1 A2^2 A3^2]

/ [A1^5 + A2^5 + A3^5

+ A1^4 A2 + A1^4 A3 + A2^4 A3

+ A1 A2^4 + A1 A3^4 + A2 A3^4

+ A1^3 A2 A3 + A1 A2*3 A3

+ A1 A2 A3^3

+ A1^2 A2^2 A3 + A1^2 A2 A3^2

+ A1 A2^2 A3^2],

where the denominator is the partition function of the problem at hand. At T=∞, one obtains E = 25(E1+E2+E3)/17, which reflects equipartitioning in the 'classical' region, and, for E1 > E2, E3, at T=0 one appropriately obtains E = 5 E1.

Incidentally, physically the relevant quantity to be considered is not the energy E, but the Helmhotlz free energy F = E - TS. For this, one needs to calculate the entropy S. This calculation is carried out along the lines of the calculation of E described above.

________

[1] J. Riordan, Introduction to Combinatorial Analysis (Dover, New York, 2002).

[2] I write down the zeros to indicate empty energy levels.

https://en.wikipedia.org/wiki/Partition_%28number_theory%29