Consider two Gaussian random variables $X\sim N(0,\sigma_X^2)$ and $Y\sim N(0,\sigma_Y^2)$ with known $E[XY] = \sigma_{XY}$. The question is how to find $E[X^2\mid Y]$ assuming $X$ and $Y$ are jointly Gaussian.
My approach: I thought of the following steps:
Find the pdf $f_Y(y)$;
Find the joint pdf $f_{X,Y}(x,y)$;
Now, we have the conditional pdf $f_{X\mid Y}(x\mid y) = f_{X,Y}(x,y)/f_Y(y)$;
$E[X^2\mid Y=y] = \int_{-\infty}^{\infty}x^2f_{X\mid Y}(x\mid y)\ dx$.
Is my approach correct? If yes, are there any other more elegant solutions?
Your approach is right if the distribution of $X$ and $Y$ is jointly Gaussian, rather than only each by itself being Gaussian. For example, suppose one had $Y = \begin{cases} +X & \text{if } |X|>c, \\ -X & \text{if } |X|<c. \end{cases}\quad$ Then it can be shown that $X$ and $Y$ are both Gaussian if $X$ is Gaussian, but the pair $(X,Y)$ is not Gaussian, since the probability that it is $0$ is positive.
Also, I would write $f_{X,Y}$ rather than $f_{XY}$ so as not to confuse it with the density of the product $XY.$
You don't need step 1.
The density is \begin{align} & c\cdot\exp \left( \frac{-1}{2(1-\rho^2)} \left( \left( \frac x {\sigma_X}\right)^2 + \left( \frac y {\sigma_Y} \right)^2 - 2\rho \left( \frac x {\sigma_X} \right) \left( \frac y {\sigma_Y} \right) \right) \right) \\[6pt] & \qquad \text{where } \rho = \frac{\sigma_{X,Y}}{\sigma_X \sigma_Y} = \operatorname{cor} (X,Y). \\[4pt] & \qquad \text{and } c = \tfrac 1 {2\pi\sqrt{\sigma_X^2\sigma_Y^2- \sigma_{X,Y}^2}} \end{align} To find the conditional density of $X$ given $Y$ we view the quadratic function of $x$ and $y$ just as a function of $x$ and complete the square: \begin{align} & \left( \frac x {\sigma_X}\right)^2 + \left( \frac y {\sigma_Y} \right)^2 - 2\rho \left( \frac x {\sigma_X} \right) \left( \frac y {\sigma_Y} \right) \\[8pt] = {} & \left[ \left( \frac x {\sigma_X}\right)^2 - 2\rho \left( \frac x {\sigma_X} \right) \left( \frac y {\sigma_Y} \right) \right] + \left( \frac y {\sigma_Y} \right)^2 \\[8pt] = {} & \left[ \left( \frac x {\sigma_X}\right)^2 - 2\rho \left( \frac x {\sigma_X} \right) \left( \frac y {\sigma_Y} \right) + \rho^2 \left( \frac y {\sigma_Y} \right)^2 \right] + \left( \frac y {\sigma_Y} \right)^2 - \rho^2 \left( \frac y {\sigma_Y} \right)^2 \\[8pt] = {} & \left[ \frac x {\sigma_X} - \rho\cdot\frac y {\sigma_Y} \right]^2 + {} \underbrace{ (1-\rho^2) \left( \frac y {\sigma_Y} \right)^2}_\text{No “$x$” appears here.} \end{align} Things not depending on $x$ are in this context constants, so we have \begin{align} f_{X\,\mid\,Y\,=\,y} (x) & = \text{constant} \times\exp\left( -\frac 1 {2(1-\rho^2)} \left[ \frac x {\sigma_X} - \rho\cdot\frac y {\sigma_Y} \right]^2 \right) \\[8pt] & = \text{constant} \times \exp\left( -\frac 1 {2\sigma_X^2(1-\rho^2)} \left[ x - \frac{\sigma_{X,Y}}{\sigma_Y^2}\cdot y \right]^2 \right) \end{align} This is a Gaussian density with variance $\sigma_X^2 (1-\rho)^2$ and expectation $\dfrac{\sigma_{X,Y}}{\sigma_Y^2} \cdot y.$
The expected value of the square of a random variable is the sum of its variance and the square of its expected value, thus it is $$ \sigma_X^2(1-\rho^2) + \left( \frac{\sigma_{X,Y}}{\sigma_Y^2} \cdot y\right)^2. $$
This can be viewed as a weighted average: $$ (1-\rho^2) \sigma_X^2 + \rho^2\left( \frac{\sigma_X}{\sigma_Y}\cdot y \right)^2. $$