Martingale Approach
As you noted, you need to solve
\begin{eqnarray}
F(0) & = & e^{-r T} \mathbb{E} \left[ \left( X_T - K \right)^2 \right]\\
& = & e^{-r T} \left( \mathbb{E} \left[ X_T^2 \right] - 2 K \mathbb{E} \left[ X_T \right] + K^2 \right)
\end{eqnarray}
Let $Y_t = X_t^2$. Then, by applying the Itō formula, we get
\begin{eqnarray}
\mathrm{d}Y_t & = & 2 X_t \mathrm{d}X_t + \mathrm{d} \langle X \rangle_t\\
& = & \left( 2 r + \sigma^2 \right) Y_t \mathrm{d}t + 2 \sigma Y_t \mathrm{d}W_t.
\end{eqnarray}
It follows that
\begin{eqnarray}
\mathbb{E} \left[ X_T \right] & = & X_0 e^{r T},\\
\mathbb{E} \left[ X_T^2 \right] & = & X_0^2 e^{\left( 2 r + \sigma^2 \right) T}.
\end{eqnarray}
Consequently,
\begin{equation}
F(0) = e^{-r T} \left( X_0^2 e^{\left( 2 r + \sigma^2 \right) T} - 2 K X_0 e^{r T} + K^2 \right).
\end{equation}
PDE Approach
Alternatively, you can apply a change of variables to the PDE. Define
\begin{equation}
\tau = (T - t), \quad \xi = r - \frac{1}{2} \sigma^2, \quad y = \ln(x) + \xi \tau, \quad F(t, x) = e^{-r \tau} G(\tau, y).
\end{equation}
After carefully applying the chain rule, we obtain the PDE
\begin{equation}
\frac{\partial G}{\partial \tau} = \frac{1}{2} \sigma^2 \frac{\partial^2 G}{\partial y^2}
\end{equation}
with initial condition
\begin{equation}
G(0, y) = \left( e^y - K \right)^2.
\end{equation}
The fundamental solution is the heat kernel given by
\begin{equation}
\phi(\tau, y) = \frac{1}{\sqrt{2 \pi \sigma^2 \tau}} \exp \left\{ -\frac{y^2}{2 \sigma^2 \tau} \right\}.
\end{equation}
We obtain $G(\tau, y)$ through the convolution
\begin{eqnarray}
G(\tau, y) & = & \int_{-\infty}^\infty G(0, z) \phi(\tau, y - z) \mathrm{d}z\\
& = & \int_{-\infty}^\infty \left( e^{2 z} - 2 K e^z + K^2 \right) \phi(\tau, y - z) \mathrm{d}z.
\end{eqnarray}
Next note that for some $\alpha \in \mathbb{R}$,
\begin{eqnarray}
\int_{-\infty}^\infty e^{\alpha z} \phi(\tau, y - z) \mathrm{d}z & = & e^{\alpha y + \alpha^2 \sigma^2 \tau / 2}
\end{eqnarray}
Thus
\begin{equation}
G(\tau, y) = e^{2 y + 2 \sigma^2 \tau} - 2 K e^{y + \sigma^2 \tau / 2} + K^2
\end{equation}
and substituting back yields
\begin{equation}
F(0, x) = e^{-r T} \left( x^2 e^{\left( 2 r + \sigma^2 \right) T} - 2 K x e^{r T} + K^2 \right)
\end{equation}
as before.
