To keep things simple, we will approximate $\nabla^2u(\textbf{x},\dot{\textbf{x}},t)$ with some operator $L[u(\textbf{x},\dot{\textbf{x}},t)]$ that uses finite difference approximations:
$$ \nabla^2u(\textbf{x},\dot{\textbf{x}},t) \approx L[u(\textbf{x},\dot{\textbf{x}},t)] = \sum_{k=1}^{d} \frac{u(\textbf{x}-h\hat{\textbf{e}}_k,\dot{\textbf{x}},t) - 2u(\textbf{x},\dot{\textbf{x}},t) + u(\textbf{x}+h\hat{\textbf{e}}_k,\dot{\textbf{x}},t)}{h^2} $$
where $d$ is the dimension of $\textbf{x}$, $h \ll 1$, and $\hat{\textbf{e}}_k$ is the unit vector in the $k^{th}$ direction. Based on this approximation, we can then view the solution to this problem as finding a path based on the constraining "dynamics" you have listed. This can then be viewed from the perspective of Calculus of Variations, where we wish to find an optimal path. This allows us to specify a cost function to minimize, such as:
$$J[\textbf{x}] = \int_{t_1}^{t_2}\left( L[u(\textbf{x},\dot{\textbf{x}},t)] - f(\textbf{x},\dot{\textbf{x}},t)\right)^2 dt$$
where $t_1$ and $t_2$ bound the time frame you care about, and your optimal path $\textbf{x}$ will minimize the cost. To make this path more easily approximated, let's approximate $\textbf{x}(t)$ by some basis, expressed as the following:
$$ \textbf{x}_{h}(t) = \sum_{j}^{n} \textbf{a}_{j} \phi_{j}(t)$$
$$ \dot{\textbf{x}}_{h}(t) = \sum_{j}^{n} \textbf{a}_{j} \dot{\phi}_{j}(t)$$
By using this approximation, $\textbf{x}_{h}(t)$, in place of $\textbf{x}(t)$, the problem can be tackled by finding the coefficients $\textbf{a}_j$ such that the cost function is minimized. This optimization looks like it will be ugly, but this is at least one way you might go about solving it.
