I've seen the nice proof of this using spheres, but I'm looking for a way to prove it parametrically if possible. Using a cylinder $x^2+y^2=r^2$ and a plane $ax+by+cz+d=0$ I got:
$x=r\cos(\theta), y=r\sin(\theta), z=\dfrac{-ar\cos(\theta)+br\sin(\theta)+d}{c}$
But after this I'm stuck trying different projections and messing with ellipse definitions
On
W.l.o.g. we can take the cylinder to have unit radius and the cutting plane to be the $x$-$y$ plane rotated through an angle of $\theta$ about the $x$-axis. This rotation is represented by the homogeneous transformation matrix $$R = \begin{bmatrix} 1&0&0&0 \\ 0& \cos\theta & -\sin\theta &0 \\ 0& \sin\theta & \cos\theta &0 \\ 0&0&0&1 \end{bmatrix}.$$ The cutting plane is represented by the homogeneous vector $\mathbf\pi = [0:-\sin\theta:\cos\theta:0]$ (these are just the coefficients of a point-normal equation of the plane). The curve of intersection can be generated by projecting the unit circle in the $x$-$y$ plane parallel to the cylinder’s axis—the $z$-axis—onto the cutting plane. If we let $\mathbf p = [\cos t:\sin t:0:1]$ be a point on the circle and $\mathbf q=[0:0:1:0]$ be the point at infinity on the $z$-axis, the projection of $\mathbf p$ onto $\pi$ is the intersection of the line $\overline{\mathbf p\mathbf q}$ with $\mathbf\pi$: $$\mathbf p' = (\mathbf\pi^T \mathbf q)\mathbf p-(\mathbf\pi^T \mathbf p)\mathbf q = [\cos\theta \cos t : \cos\theta \sin t : \sin\theta \sin t : \cos\theta].$$ Now rotate this by $R^{-1}=R^T$ to bring it back onto the $x$-$y$ plane. This produces $[\cos\theta \cos t : \sin t : 0 : \cos\theta]$, which is the point with inhomogeneous coordinates $\left(\cos t,\sec\theta \sin t,0\right)$. If you don’t recognize this as the parametrization of an ellipse with semimajor axis length $\sec\theta$ and semiminor axis length $1$, you can eliminate $t$: $$x^2 + y^2\cos^2\theta = \cos^2t+\sin^2t = 1.$$
Without loss of generality let's assume the cylinder has radius $1$, so that it has equation $x^2 + y^2 = 1$. Assuming that the plane is not parallel to the cylinder, we can always rearrange the coordinate system so that the plane goes through the origin, or even better, make the plane go through the $x$-axis after a suitable rotation. Now the plane should have equation $z = y \tan\alpha$ (with $\alpha$ being the slope in the $yz$-plane).
Now your parametrization gives the curve $$x = \cos\theta, \quad y = \sin\theta, \quad z=\sin\theta\tan\alpha.$$ We think of a rotation (of the whole space) around the $x$-axis, through an angle of $-\alpha$ in the $yz$-plane: $$\begin{aligned} \begin{bmatrix} x \\ y \\ z \end{bmatrix} & \mapsto \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos(-\alpha) & -\sin(-\alpha) \\ 0 & \sin(-\alpha) & \cos(-\alpha) \end{bmatrix} \begin{bmatrix} x \\ y \\ z \end{bmatrix} \\ & = \begin{bmatrix} x \\ y\cos\alpha + z\sin\alpha \\ -y\sin\alpha + z\cos\alpha \end{bmatrix}. \end{aligned}$$ Under this rotation the curve becomes $$\begin{aligned} \begin{bmatrix} \cos\theta \\ \sin\theta \\ \sin\theta\tan\alpha \end{bmatrix} & \mapsto \begin{bmatrix} \cos\theta \\ \sin\theta\cos\alpha + \sin\theta\tan\alpha\sin\alpha \\ -\sin\theta\sin\alpha + \sin\theta\tan\alpha\cos\alpha \end{bmatrix} \\ & = \begin{bmatrix} \cos\theta \\ \sin\theta\cos\alpha + \sin\theta\sin^2\alpha/\cos\alpha \\ -\sin\theta\sin\alpha + \sin\theta\sin\alpha \end{bmatrix} \\ & = \begin{bmatrix} \cos\theta \\ \sin\theta/\cos\alpha \\ 0 \end{bmatrix}. \end{aligned}$$ This is $$x = \cos\theta, \quad y = \frac{1}{\cos\alpha}\sin\theta, \quad z\equiv0,$$ which is an ellipse.