I still don't understand how to get your formula for $n$. If I understand correctly, you assimilate each rod by its exclusion volume, giving it three degrees of freedom, $2$ of translation, one of rotation. Let $A=\theta$ the space occupied by each effective rod in rotation space (rotation within the excluded volume), and $V$ the total accessible position space, I get the number of microstates up to a $\theta$ independent factor:
$$
W = A^NV(V-\Omega) ...(V-(N-1)\Omega)\\
\propto A^N\prod_{k = 0}^{N-1}\left(1-n\Omega\frac{k}{N}\right)
$$
($\theta$ independent factor) In the thermodynamic limit $N\gg1$ this becomes:
$$
S =N\ln A +N\int_0^1dx\ln(1-n\Omega x) +S_0\\
= N\left(\ln A -1-\frac{1-n\Omega}{n\Omega}\ln(1-n\Omega)\right) + S_0
$$
($S_0$ a $\theta$ independent additive term)
so
$$
\frac{dS}{d\theta} = N\frac{A'}{A}+N\left(\frac{1}{n\Omega}+\frac{1}{(n\Omega)^2}\ln(1-n\Omega)\right)n\Omega'
$$
From the equation $\frac{dS}{d\theta} = 0$ you have, with some caveats, an implicit definition of $\theta$. Btw, there is no simple formula for $n$ as a function of $\theta$. More precisely, for all of this to make sense, you need to assume all along that $n\Omega\leq 1$ so letting $n_c = \frac{1}{\pi l^2}$, for $n\leq n_c$ the above resolution is for $\theta\in[0,\pi]$ and for $n\geq n_c$ you have $\theta\in[0,\theta_m]$ with the implicit definition $\Omega(\theta_m) = 1/n$.
For the actual resolution of $\theta$, you therefore need to distinguish the two cases:
$n>n_c$: there is necessarily have a zero since $\frac{dS}{d\theta} \to +\infty$ when $\theta\to 0$ and $\frac{dS}{d\theta} \to -\infty$ when $\theta\to \theta_m$. It is actually the only solution, which defines $\theta$.
$n<n_c$: there is no solution, and you always have $\frac{dS}{d\theta}>0$, so $\theta=\pi$.
On an intuitive level, you have a competition between the number of states in rotation which increases with $\Omega$ and number of states in position which decreases with $\Omega$. When $n\gg n_c$, the second effect dominates and $\theta$ is small and as $n$ decreases, $\theta$ increases. At high density, the rods are tightly packed and are therefore highly constrained. Conversely, at low density, the rods are completely free to rotate. This is why you get a phase transition, you have a sudden plateau at $n=n_c$ and a resulting discontinuous slope. This kind of plateau effect is recurrent in those self-consistent approaches, you might have seen something similar in mean-field theory.
You can look at the asymptotic behaviour of $\theta$ when $n\to n_c^+$. Let $n = n_c+\epsilon$ and $\theta=\pi-\delta$ in this regime, I found (the computation is a bit tricky, you should check it for yourself):
$$
\delta = (2\pi^2\epsilon)^{1/3}
$$
which gives you a discontinuous derivative at the transition.
Hope this helps.
