I have posted a similar question on math.stackexchange (https://math.stackexchange.com/questions/1848492/calculate-mean-of-sde), but didn't find anyone who could help.
I'm interested in the one-dimensional SDE
$$ \mathrm d u = -u^3 \mathrm d t + u^2 \mathrm d W$$
where $u$ is adapted to the Brownian motion $W$. By a Lyapunov function argument (and by local Lipschitz continuity), we can show that unique solutions exist and do not blow up in finite time with probability 1.
To get a feel for the qualitative behavior, here's a simulation of its paths. The first picture shows sample paths and the empirical mean, the second picture shows the empirical mean and the empirical third moment.
According to the Ito formula (see the posted linked above), $\mathbb E u^3(t) = \mathbb E u^3(0)$. This can be seen in the simulations (although the number of samples is not high enough: The third moment stays roughly constant for a short time but then declines because the rare events needed to have constant third moment do not happen in the Monte Carlo simulation. For short time this looks fine, though).
Now my problem is the mean of the differential equation:
$$ \begin{align*} u(t) &= u(0) + \int_0^t u^3(s)\mathrm d s + \int_0^t u^2(s) \mathrm d W_s\\ \mathbb E u(t) &= \mathbb E u(0) + \mathbb E\int_0^t u^3(s)\mathrm d s\\ &= \mathbb E u(0) + t\cdot \mathbb E u^3(0) \end{align*} $$
which obviously cannot be correct (see the empirical mean which gives a completely different example). There are a few fishy steps in that calculation, for example the unjustified swapping expectation and integral, but I can't see how this could give such a blatantly wrong result.${}$