Residuals
Given the data $\lbrace \mathbf{x_1}, \mathbf{x_2}, ... , \mathbf{x_n} \rbrace$ you can compute for each model the $n-1$ residuals:
$$\mathbf{\epsilon_{i+1}} = f(\mathbf{x_i}) - \mathbf{x_{i+1}}$$
and the likelihood is computed according to $\mathbf{\epsilon} \sim N(\boldsymbol{\mu},\boldsymbol{\Sigma})$
Simpler expression to use in least sum of squares
Possibly you have multiple models, or you want to optimize some parameter.
When you set $\mathbf{y_i} = \mathbf{x_{i+1}}$ then your problem can be expressed more simple like:
$$\mathbf{y_i} = f(\mathbf{x_{i}})+\mathbf{\epsilon_i}$$
This you can solve with any of the multitude of least sum of squares techniques.
Sidenote
Your model does assume that the randomness in the observations is not in a measurement error (or that those errors are negligible) but only in the generation process of the next $\mathbf{x_{i+1}}$ from the past $\mathbf{x_{i}}$.
If your process is related to true variables $\mathbf{x_{i+1}}$ and $\mathbf{x_{i}}$
$$ \mathbf{x_{i+1}} = f(\mathbf{x_{i}}) + \mathbf{\epsilon_{i+1}}$$
but your observed variables $\mathbf{y_{i}}$ have errors in relation to those true (underlying) variables $\mathbf{y_{i}} = \mathbf{x_{i}} + \mathbf{\epsilon^\prime_{i}}$ then your model to fit (which relates the $\mathbf{y_{i}}$ and not the $\mathbf{x_{i}}$) should be
$$ \mathbf{y_{i+1}} + \mathbf{\epsilon^\prime_{i+1}} = f(\mathbf{y_{i}}+\mathbf{\epsilon^\prime_{i}}) + \mathbf{\epsilon_{i+1}}$$
since you say that $f(\cdot)$ is a matrix you can solve this (it would be more difficult when $f(\cdot)$ is a non-linear function), and you get a model like before $\mathbf{y_i} = f(\mathbf{x_{i}})+\mathbf{\epsilon_i}$ but the succeeding errors $\mathbf{\epsilon_{i}}$ and $\mathbf{\epsilon_{i+1}}$ are correlated
