your final answer isn't right

You're right. That should be ε, not E. I did the extra few steps to substitute α = E/(Nε) back in, and solve for E, to recover DanielFilan's (corrected) result:

E = Nε / (exp(ε/T) + 1)

I used S = log[N choose M], where M is the number of excited particles (so M = αN). Then I used Stirling's approximation as you suggested, and differentiated with respect to α.