You're right on track but don't have enough equations. The equations which give you the solution are:
\begin{align}
k&=\frac{Z_1 Z_2 e^2}{4 \pi \varepsilon_0} \tag{1}\\
E&=\frac{1}{2}m v_{min}^2+\frac{k}{r_{min}} \tag{2}\\
\frac{1}{2}m v_{min}^2&=E \frac{b^2}{r_{min}^2} \tag{3}\\
b&=\frac{k}{2 E} \cot\left(\frac{\theta}{2}\right)\tag{4}
\end{align}
Equation 1 is just shorthand. Equation 2 is your second equation, but with $\frac{1}{2} mv_0^2$ replaced with $E$. Equation 3 is your fourth equation ($v_{min}=\frac{v_0 b}{r_{min}}$), but squared, multiplied by $\frac{1}{2}m$, and with $\frac{1}{2} mv_0^2$ replaced by $E$. Equation 4, finally, is the important result of Rutherford scattering relating the impact parameter to the scattering angle.
Because I'm lazy I threw it into Mathematica, but if you make these substitutions (most importantly: Getting $v_{min}$, $v_0$, and $b$ OUT of your quadratic equation), it should be easy to do by hand.
The mathematica code
FullSimplify[
Solve[{eE == 1/2 m vmin^2 + k/rmin, 1/2 m vmin^2 == eE b^2/rmin^2,
b == k/(2 eE) Cot[\[Theta]/2]}, {rmin, vmin, b}],
Assumptions -> {eE > 0, m > 0}]
spits out, as one solution,
$$\frac{k \left(\sqrt{\csc ^2\left(\frac{\theta }{2}\right)}+1\right)}{2 \text{eE}}$$
which is exactly what you need. ("eE" is a single variable used in the mathematica code, and is not a product. This is because mathematica hates single capital letters. It represents the energy of the system, $\frac{1}{2} m v_0^2$.)
In the Wikipedia article about rutherford scattering the derivation of the scattering cross section
$$ \frac{d\sigma}{d\Omega} =\left(\frac{ Z_1 Z_2 e^2}{8\pi\epsilon_0 m v_0^2}\right)^2 \csc^4{\left(\frac{\Theta}{2}\right)}$$
is given. Let's rewrite that in your notation: $Z_1 = Z$, $Z_2 = 4$, $k = \frac{1}{4\pi\varepsilon_0}$ and $KE = \frac{1}{2}m v^2$:
$$ \frac{d\sigma}{d\Omega} =\left(\frac{ Z e^2 k}{KE}\right)^2 \csc^4{\left(\frac{\Theta}{2}\right)} = \frac{Z^2 k^2 e^4}{ KE^2\sin^4(\theta/2)}$$
The relation between the cross section and the number of detected particles is
$$\frac{N(\theta)}{N_i} = \frac{nL}{4r^2} \frac{d\sigma}{d\Omega}$$
$nL$ gives the number of targets to scatter from and $1/(4r^2)$ is because at a bigger distance the intensity goes down by this factor (I'm sure about the $r^2$ but not exactly about the factor. I thought it would be $2\pi$ so maybe I did an error there.)
I hope this helps.
Best Answer
$N(\theta)/N_i$ is a probability density not a probability. To get a probability you have to integrate it over $\theta$ or maybe multiply it by $d\theta$. That makes the numbers reasonable again.
Usually the difference is clear because of dimensions, if you're dealing with distributions in things like time or position or energy. But angles are dimensionless that doesn't help here.