We have an expression for $\Omega$ (occupation of each macrostate) in terms of $n_i$ (occupation numbers) . We want to find the $n_i$ which maximises $\Omega$.
We now that
$$ln[\Omega]=ln[N!]-\sum_i(ln[n_i]!)$$
Using stirlings approximation we can write this as
$$=(Nln[N]-N)-\sum_i(n_iln[n_i]-n_i)$$
We are subject to two constraints 1. Total number of particles remain constant
$$N=constant: \sum_idn_i=0$$
- Total internal energy remains constant
$$U=constant: \sum_i\epsilon_i dn_i=0$$
We now use lagrange's method of undetermined multipliers.
$$d(ln[\Omega]=-d(\sum_i(n_iln[n_i]-n_i))=0$$
$$=-\sum_i dn_i(ln[n_i]+\frac{n_i}{n_i}-1)=0$$
This implies that
$$\sum_iln[n_i]dn_i=0$$
So, we have three things which must be zero simultaneously. We can multiple zero by any factor and still have zero so,
$$\sum_iln[n_i]dn_i=0, \alpha\sum_i\epsilon_i dn_i=0, \beta\sum_idn_i=0 $$
We can combine these expressions to write our overall maximisation problem as:
$$\sum_i(ln[n_i]+\alpha+\beta\epsilon_i)dn_i=0 \forall \alpha,\beta$$
Now the part I raise issue with is: apparently the only way for this to be solved is that every term in the sum has to be zero. So
$$ln[n_i]+\alpha+\beta\epsilon_i=0 \forall i,$$ but specific $\alpha$ and $\beta$.
However, surely, we are doing a sum here, there is no requirement that this be zero at all points. For instance, the sum may be zero, but part of the sum might be positive, but this is canceled by a negative part.
Anyway, just to finish the derivation, this rearranges to
$$ln[n_i]=-\alpha-\beta\epsilon_i$$
Solving for $n_i$
$$n_i=e^{-\alpha-\beta\epsilon_i}=Ae^{-\beta\epsilon_i}$$
Giving
$$N=\sum_in_i=A\sum_i^{\beta\epsilon_i}$$