You are dealing with multiple sources of errors. Namely, the ODE solver ode45 does not give an exact solution for $P(t)$ but a numerical approximation for it and the method you use to simulate the state is forward Euler, which is only a first order method. Since your system is unstable any perturbation might get amplified a lot. Also keep in mind the cost function that you are trying to minimize, because it might not be cost effective for finite horizon LQR to drive the state close to zero.
Fortunately one can improve the accuracy of many of the numerical results. Namely, there is an analytical solution for $P(t)$, and $x(t)$ can also be solved with ode45. The analytical solution for $P(t)$ can be found by using $P(t) = \bar{P} + V^{-1}(t)$, with $\bar{P}$ the stationary solution of $P(t)$ so by solving the related CARE. The dynamics of $V(t)$ can shown to be
$$
\dot{V} = V\,\mathcal{A}^\top + \mathcal{A}\,V - B\,R^{-1} B^\top,
$$
with $\mathcal{A} = A - B\,R^{-1} B^\top \bar{P}$. By using vectorization and the Kronecker product that dynamics can also be written as
$$
\dot{z} = M\,z,
$$
with $M = I \otimes \mathcal{A} + \mathcal{A} \otimes I$ and $z(t) = \text{vec}(V(t)) - M^{-1}\,\text{vec}(B\,R^{-1} B^\top)$. The solution for $z$ is given by
$$
z(t) = e^{M\,(t - \tau)}\,z(\tau).
$$
So $P(t)$, given $P(T)$, can be obtained using
\begin{align}
V(T) &= \left(P(T) - \bar{P}\right)^{-1}, \\
z(T) &= \text{vec}(V(T)) - M^{-1}\,\text{vec}(B\,R^{-1} B^\top), \\
V(t) &= \text{vec}^{-1}\!\left(e^{M\,(t - T)}\,z(T) + M^{-1}\,\text{vec}(B\,R^{-1} B^\top)\right), \\
P(t) &= \bar{P} + V^{-1}(t).
\end{align}
I also did some simulations using your system/script while changing only your numerical solution for $P(t)$ with the analytical solution, the integration method for the dynamics of the state from forward Euler to ode45 or both. From those results it does seem that using the analytical solution $P(t)$ seemed to have to biggest impact.
Consider the following optimal control problem
$$
J_T(x_0)=\min_{u}\sum_{i=0}^{T-1}x_i^TQx_i+u_i^TRu_i
$$
where $x_{k+1}=Ax_k+Bu_k$. From dynamic programming, we have that
$$
J_T(x_0)=x_0^TP_T(0)x_0
$$
where $P_T(k)$ solves the Difference Riccati Equation (DRE)
$$
P_T(k-1)=\mathcal{R}(P_T(k)),\ P_T(T)=0.
$$
In that case, the optimal control law is given by $u_k=K_kx_k$ where $K_k=-(R+B^TP_T(k+1)B)^{-1}B^TP_T(k+1)A$.
On the other, we can define the cost
$$
J_\infty(x_0)=\min_{u}\sum_{i=0}^{
\infty}x_i^TQx_i+u_i^TRu_i.
$$
and we know that this cost is equal to
$$
J_\infty(x_0)=x_0^TP_\infty x_0
$$
where $P_\infty$ is the stabilizing solution of the Algebraic Riccati Equation (ARE)
$$
P_\infty=\mathcal{R}(P_\infty).
$$
By stabilizing solution, it is meant here that $A+BK$ is Schur stable where
$K=-(R+B^TP_\infty B)^{-1}B^TP_\infty A$.
The following still contains gap and unclear points which I will fix when I understand how to do it correctly. In the meantime that can just serve as an indicative answer. We show that $P_N(0)\to P_\infty$ as $T\to\infty$.
First of all, it is important to note that $P_{T+1}(0)\succeq P_T(0)$ because the integrand is nonnegative. Assumming now that the pair $(A,B)$ is stabilizable, there exists a $K$ such that $A+BK$ is Schur stable.
This means that we have
$$
J_T(x_0)\le\min_{u}\sum_{i=0}^{T-1}x_i^T(Q+K^TRK)x_i=x_0^T\left(\sum_{i=0}^{T-1}((A+BK)^i)^T(Q+K^TRK)A+BK)^i\right)
$$
and we have that the sum is converging as $T\to\infty$ to, say, $Q$. Therefore, this implies that $P_T(0)\preceq Q$. As a result, the sequence $\{P_T(0)\}_T$ is nondcreasing and upper bounded, which means that it is converging to a limit $P$, which is assumed to be positive definite.
Since it is the solution of a sequence, then it must be converging to a stationary solution of the DRE. In fact, since it is positive definite, it is the unique positive definite solution to the ARE and we have that $P=P_\infty$ (this needs to be better supported).
Best Answer
We have that
$$ \begin{array}{rcl} C(K)&=&\mathbb{E}\left[\sum_tx_t^TQx_T+u_t^TRu_t\right]\\ &=&\sum_t\mathbb{E}\left[x_t^T(Q+K^TRK)x_t\right]\\ &=&\mathrm{trace}\left[(Q+K^TRK)\mathbb{E}\left[\sum_t x_tx_t^T\right]\right]\\ &=&\mathrm{trace}((Q+K^TRK)P). \end{array} $$ where we have set $P:=\mathbb{E}\left[\sum_tx_tx_t^T\right]$.
If we let $S_t:=\mathbb{E}[x_tx_t^T|S_{t-1}]$, then we have that
$$ \begin{array}{rcl} S_{t+1}&=&\mathbb{E}[x_{t+1}x_{t+1}^T|S_{t}]$\\ &=&(A+BK)S_t(A+BK)^T \end{array} $$ and $S_0=\mathbb{E}[x_0x_0^T]=I$. So, we have that
$$ S_t=(A+BK)^t((A+BK)^T)^t. $$
It can be shown that $P$ exists if and only if $A+BK$ is Schur stable, which is the case here. This yields
$$ P=\sum_{t=0}^\infty(A+BK)^t((A+BK)^T)^t. $$ We can play around this expression as $$ \begin{array}{rcl} P &=& I+\sum_{t=1}^\infty(A+BK)^t((A+BK)^T)^t\\ &=& I+(A+BK)\left[\sum_{t=1}^\infty(A+BK)^{t-1}((A+BK)^T)^{t-1}\right](A+BK)^T\\ &=& I+(A+BK)\left[\sum_{t=0}^\infty(A+BK)^{t}((A+BK)^T)^{t}\right](A+BK)^T\\ P &=& I+(A+BK)P(A+BK)^T \end{array} $$ to get the desired expression $$ (A+BK)P(A+BK)^T-P+I=0. $$