Consider a simple sde $$\mathrm{d}X_t = \mu \mathrm{d}t + \sigma \mathrm{d}W_t.$$ Assuming $X_0 = x$, $$\mathbb{E}\left[X^2_t\right] = \left(x +\mu t\right)^2 + \sigma^2 t.$$ I have been seeking to approximate this curve by sampling the SDE. I employ an Euler-Maruyama method to generate M discretized paths. To estimate $\mathbb{E}\left[X^2_t\right]$, I fix a time point, then average the square value of each of the paths at the fixed time point.
The resulting curve does not agree with the analytical solution. It is exponentially greater than the correct curve.
Is there something I'm missing about my method? If, for example, I wanted to approximate the solution to a pde using Feynman-Kac, I would essentially do the same thing for $u(t,x) = \mathbb{E}\left[f\left(X_t\right)\middle|X_0 = x\right]$. For a fixed $x$ in this case, I can recover the correct $u(t,x)$ using the exact same method. This leads me to believe I should be able to do the same thing when $f(x) = x^2$.
Is there something strange about the second moment that I'm missing?