Proof that any plane perpendicular to x,y plane intersection with the bivariate normal distribution has the shape of a normal distribution. The bivariate normal distribution is given by the equation:
$$f(x,y)=\frac{\exp\left(-\frac{1}{2(1-\rho)^2}\left[\left(\frac{x-\mu_1}{\sigma_1}\right)^2-2\rho\left(\frac{x-\mu_1}{\sigma_1}\right)\left(\frac{y-\mu_2}{\sigma_2}\right)+\left(\frac{y-\mu_2}{\sigma_2}\right)^2\right]\right)}{2\pi\sigma_1\sigma_2\sqrt{1-\rho^2}}$$
for $-\infty<x<\infty$ and $-\infty<y<\infty$, where $\sigma_1>0$, $\sigma_2>0$, and $-1<p<1$.
The normal distribution is given by the equation:
$$f(x)=\frac{e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2}}{\sigma\sqrt{2\pi}}$$
for $-\infty<x<\infty$ where $\sigma>0$.



I am trying to show that if I take any plane that is perpendicular to the $x,y$ plane, it's intersection with the bivariate normal distribution  is equivlent to $c\cdot f(x)$, where $f(x)$ is the normal distribution and $c$ is a real constant.
This can also be written as showing the following statement to be true.
$$\frac{\exp\left(-\frac{1}{2(1-\rho)^2}\left[\left(\frac{x-\mu_1}{\sigma_1}\right)^2-2\rho\left(\frac{x-\mu_1}{\sigma_1}\right)\left(\frac{(mx+b)-\mu_2}{\sigma_2}\right)+\left(\frac{(mx+b)-\mu_2}{\sigma_2}\right)^2\right]\right)}{2\pi\sigma_1\sigma_2\sqrt{1-\rho^2}}=c\cdot \frac{e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2}}{\sigma\sqrt{2\pi}}$$
or
$$f(x,mx+b)=c\cdot f(x)$$
Where $m$ and $b$ are any real number.
This does leave out all the planes parallel to the $x$ axis, although those will be trivial to show separately, after showing that $f(y,x)$ is also a bivariate normal distribution.
I tried simplifying this but have made no progress so far, It's quite likely that I am going about this in completely the wrong way.

As you can see from another answer, it is possible to carry through
your initial idea to get a proof.
The one detail that I questioned was how we establish that the
coefficient of $x^2$ in
$$
\frac{1}{(1-\rho)^2}\left[\left(\frac{x-\mu_1}{\sigma_1}\right)^2-2\rho\left(\frac{x-\mu_1}{\sigma_1}\right)\left(\frac{mx+b-\mu_2}{\sigma_2}\right)+\left(\frac{mx+b-\mu_2}{\sigma_2}\right)^2\right]
$$
is positive (using $m$ and $b$ as defined in the question). In fact, that coefficient is
$$
\alpha = 
\frac{1}{(1-\rho)^2} \left(\frac{1}{\sigma_1^2}-2\rho\frac{m}{\sigma_1\sigma_2}+\frac{m^2}{\sigma_2^2}\right) .
$$
But
$$
\frac{1}{\sigma_1^2}-2\rho\frac{m}{\sigma_1\sigma_2}+\frac{m^2}{\sigma_2^2}
= \left(\frac{1}{\sigma_1}-\rho\frac{m}{\sigma_2}\right)^2 
+  (1 - \rho^2)\frac{m^2}{\sigma_2^2},
$$
and since $-1<\rho<1$ it follows that $1-\rho^2>0$ and
also that $\alpha > 0.$
Note that for any $\rho$ such that $\lvert\rho\rvert>1$ 
it is possible to choose the other parameters so that $\alpha<0,$
and of course the entire polynomial is undefined for $\lvert\rho\rvert=1,$
so the fact that $-1<\rho<1$ is a necessary condition for this
proof to go through.

An alternative is a suitable substitution of variables for both $x$ and $y$ that transform the plane so that the distribution over the new variables is standard normal with covariance zero over the transformed variables and such that the equation of the line in the transformed variables has the form $y'=k.$
Such a substitution certainly exists, but it might involve
translation (to eliminate the means), rotation (to eliminate the covariance), scaling by unequal factors along the transformed
axes (to eliminate the variances), and then a final rotation to
make the perpendicular plane in the question intersect the 
transformed coordinate plane in a "horizontal" line.
This method is conceptually simple, and avoids the concern about signs
that I had with the other approach, but you may find that actually
working through the math of the other approach is easier for you.
