The argument below is paraphrased from Zangwill, Section 10.2.2. Set the origin of our coordinates at the spot we're trying to find the field, and let the axis of the solenoid be the $z$-axis. Consider a horizontal "slice" of the solenoid of width dz at a height $z$ above the xy-plane. By the Biot-Savart Law, the magnetic field due to this slice is just that of a current loop:
$$
d\vec{B} = - \frac{\mu_0 K dz}{4 \pi} \oint \frac{d\vec{\ell} \times \hat{r}}{r^2} = - \frac{\mu_0 K dz}{4 \pi} \oint \frac{d\vec{\ell} \times \vec{r}}{r^3}.
$$
(The minus sign is there because $\vec{r}$ in this formula denotes the source point, not the field point.) The slice can be thought of as a parametric curve $\vec{\ell}(\lambda)$, and it can be related to $\vec{r}$ by $\vec{r} = \vec{\ell} + z \hat{z}$, where $\vec{\ell}$ is in the xy-plane. Since $\vec{\ell}$ is perpendicular to $\hat{z}$ by construction, we have $r = \sqrt{\ell^2 + z^2}$, and so we can write the full integral as
\begin{equation}
d\vec{B} = - \frac{\mu_0 K dz}{4 \pi} \oint \frac{d\vec{\ell} \times (\vec{\ell} + z \hat{z})}{(\ell^2 + z^2)^{3/2}} = - \frac{\mu_0 K dz}{4 \pi} \left[ \oint \frac{d\vec{\ell} \times \vec{\ell} }{(\ell^2 + z^2)^{3/2}} + z \oint \frac{d\vec{\ell} \times \hat{z} }{(\ell^2 + z^2)^{3/2}} \right ]
\end{equation}
When we integrate from $z = -\infty$ to $\infty$, the second term will vanish (because it is odd with respect to $z$, while the first term evaluates to $2/\ell^2$. Thus, we now have
$$
\vec{B} = -\frac{\mu_0 K}{2 \pi} \oint \frac{d\vec{\ell} \times \vec{\ell} }{\ell^2}.
$$
Writing $\vec{\ell}$ in cylindrical coordinates, we have $\vec{\ell} = s \hat{s}$ and $d\vec{\ell} = ds \hat{s} + s \, d\phi \hat{\phi}$; thus, $-d\vec{\ell} \times \vec{\ell} = -s^2 \, d\phi \hat{z}$, and we have
$$
\vec{B} = \frac{\mu_0 K}{2 \pi} \hat{z} \oint \, d\phi.
$$
The integral is now the net solid angle traversed when we go around the loop. If the curve $\vec{\ell}(\lambda)$ encloses the origin (i.e., we are inside the solenoid), this will be $2\pi$; otherwise, it will be zero. Thus,
$$
\vec{B} = \begin{cases} \mu_0 K \hat{z} & \text{inside} \\
0 & \text{outside}.\end{cases}
$$
Note that this derivation did not assume any particular shape for the cross-section of the solenoid (i.e., the shape of the curve $\vec{\ell}$.)
To preemptively address any concerns about infinite solenoids being unrealistic: how would this derivation be modified if we had a non-infinite (but long) solenoid? Return to the step where we had written $d\vec{B}$ in terms of $dz$ and two integrals. If we assume that the solenoid now stretches from $z_1$ to $z_2$ instead, the expression becomes:
$$
\vec{B} = - \frac{\mu_0 K}{4 \pi} \left[ \oint \frac{\sigma(z) \, d\vec{\ell} \times \vec{\ell} }{\ell^2 \sqrt{(\ell/z)^2 + 1}} - \oint \frac{d\vec{\ell} \times \hat{z} }{|z| \sqrt{(\ell/z)^2 + 1}} \right ]_{z = z_1}^{z_2}.
$$
where $\sigma(z)$ is the sign of $z$. While we can't evaluate these integrals in as much generality, we still retain some important features. The first term will still point in the z-direction regardless of $z_1$ and $z_2$, and the second term will always point in the xy-plane. Moreover, the second term still vanishes so long as $z_1 = - z_2$. Thus, at the midpoint of a long solenoid, the external field still points along the axis of the solenoid (as it must by symmetry.)
If we wanted, we could expand these integrals in a power series in $\ell/z_1$ and $\ell/z_2$, assuming this quantity is small; this is saying that the perpendicular distance from the field point to the solenoid is much smaller than the distance to its ends. In this limit, we would recover the infinite-solenoid result plus corrections due to the finite length of the solenoid. The leading-order correction to the horizontal field would be $\mathcal{O}(\ell/z)$, while the leading-order correction to the vertical field will be $\mathcal{O}(\ell/z)^2$.
Best Answer
A clarification first. If we are talking about a solenoid of infinite length, then indeed the magnetic field outside is equal to zero. But for a real world solenoid its just very weak compared to the inside of the solenoid.
Let's talk about the first case first. Consider a cross-section that cuts the solenoid with a plane parallel to it. So, if you visualize it, you will see (for example) in the upper side of the cross-section the current popping out and in the lower side the current in the opposite way. Because of the infinite length, each side(the upper and lower) will produce magnetic fields that are opposite to each other, will only have horizontal direction(as each wire's vertical component of its produced magnetic field cancels out with the vertical component of the produced magnetic field of the current in the neighboring loops) and will not be dependent on the distance from the center of the solenoid(you can show that with the Biot-Savart law). The last fact is a consequence of the infinite length(infinities tend to give such "nice" results). Keeping these things in mind, the produced magnetic field from the upper and lower sides of the cross section will add up on the inside, thus giving it double the intensity while they completely cancel out(due to opposite directions) on the outside of the solenoid(due to the fact that both fields are equal in size at each point in space).
Now, for a finite-length solenoid, the magnetic fields of the upper and lower part of the cross-section will be dependent on the distance from the center of the solenoid and thus will not completely cancel out. Now, the field inside will be close but not exactly equal to double the amount of each field's strength(if the radius of the solenoid is small. If not, then the field will again be the algebraic sum of the two magnetic fields, but now the fields will not have nearly identical values at each point due to the fact that the distance of the two fields from the center will have a big difference). With the same logic, the field outside is close to zero but not zero. Also, at some point in space outside the solenoid, there will also be vertical components of the magnetic field due to the fact that near and at the edges of the solenoid, the vertical components of the magnetic field will not cancel out because there are no loops that carry current(and thus produce magnetic field) after the edged.
NOTE: In my opinion, the argument that the magnetic field outside the solenoid is weak because the magnetic field lines spread out to a very large space(out to infinity) is wrong. The reason for this is because magnetic field lines are not like lines that are used to describe the motion of water because water has limited quantities. For the description of water, you can use this as an argument for water density, because if the limited amount of water spreads out(thus the lines will be spread out) it means that the density will be low. But here, the magnetic field lines describe a field which depends a great deal of the current, so in that respect there is no shortage of it. For example, if we measure the distance from the upper side of the cross-section(z=0 at the upper side), then at z=a and z=-a (where a