Multivariate multiple non linear regression I know there are already two questions on this topic, but neither has an answer.  
I have a set of $N$ experiments. For each experiment, denoted by a vector of predictors $\mathbf{x}$, I measure $m$ values $y(f_1)\dots y(f_m)$ of a quantity $y$ at increasing frequency values. In other words, the result of each experiment is not a scalar but a function of frequency, measured at $m$ different frequencies. I thus have a total of $N\times m $ measured values of $y$.
I have a model for the results of this experiment, which is given by a computer code $g(f, \mathbf{x}, \boldsymbol{\theta})$, where the  $\boldsymbol{\theta}$ are calibration parameters. The setting is very similar to this question, with the added complication that now the output of the code is multivariate (a vector, instead than a scalar).
I would like to calibrate my code. I think I could use multivariate multiple nonlinear   least squares . However, I'm only familiar with multiple nonlinear least squares (single output). How do you calibrate on all outputs at once? Also, given a specific experiment, I'm not sure that errors at different frequencies are independent. As a matter of fact, based on past investigations I expect that (before calibration) the code results will match the trend in frequency of the measurements, but not the absolute values. What should I do in this case? 
 A: Let me change your notation just a little to make answering easier.  $y_{im}$ is the value of the outcome for experiment $i$ at frequency $m$, $i=1,\ldots,N$ and $m=1,\ldots,M$.  You have a function, $h(x,\theta)$, which outputs a $M$-dimensional vector of predicted $y$, one for each frequency, given a particular value for predictor variables, $x$, and parameter value $\theta$.  Connecting to your notation, $h(x,\theta)=g(f,x,\theta)$ (I think).  It's going to be convenient to think of a whole experiment's observations at once, so let $Y_i$ be the $M$-vector of all of experiment $i$'s $y$ for the various frequencies stacked up.  
The data-generating process is then hypothesized to be:
\begin{align}
Y_i &= h(x_i,\theta)+\epsilon_i
\end{align}
Measurement error (or whatever the source of error is) for experiment $i$ is described by the random variables $\epsilon_i$.  All of the objects in the above equation are $M$-vectors.  You say that different frequencies' data are likely to be correlated.  I assume it is also the case that different frequencies' data could have different variances.  If we are safe to assume that these variances and covariances are the same across experiments (and that each experiment is independent of each other experiment), then we can assume that $V(\epsilon_i)=\Sigma$, where $\Sigma$ is an $M \times M$ variance-covariance matrix, and $V(\epsilon)=I \otimes \Sigma$, where $\epsilon$ is the "stacked-up" vector of errors ($NM$ by 1) and $\otimes$ is Kronecker product.
The above set-up is called non-linear seemingly unrelated regression.  The standard reference in Economics is Gallant, AR (1975) Seemingly unrelated nonlinear regressions. Journal of Econometrics, 3(1): 35-50.
The way you estimate the parameters (if you know $\Sigma$) is by solving:
\begin{align}
min_\theta \sum_{i=1}^N \left( Y_i-h(X_i,\theta)\right)' \Sigma^{-1} \left( Y_i-h(X_i,\theta)\right)
\end{align}
The variance of the (consistent and asymptotically normal) estimator, $\hat{\theta}$, so-defined is:
\begin{align}
V(\hat{\theta}) &= \left( \sum_{i=1}^N \frac{\partial h}{\partial \theta}' \Sigma  \frac{\partial h}{\partial \theta} \right)^{-1}
\end{align}
Normally, you don't know what $\Sigma$ is, however, so you have to estimate it.  The way you do this is by picking a trial value of $\Sigma$ (like the identity matrix).  Then you do the minimization above to get a first try at an estimator of $\theta$, call it $\tilde{\theta}$.  Then you "estimate" the $\epsilon$ by the residuals, $e_i=Y_i-h(X_i,\tilde{\theta})$ and calculate their sample variance matrix:
\begin{align}
\hat{\Sigma} &= \frac{1}{N} \sum_{i=1}^N e_i e_i'
\end{align}
Finally, you use $\hat{\Sigma}$ in place of $\Sigma$ above in the minimization and the variance-covariance matrix.  This procedure is called two-stage, non-linear seemingly unrelated regression (or words to that effect).
In R, I think the function you want is called nlsystemfit which has documentation here.
