I have a model which takes 3 independent variables, 4 parameters and returns 2 values.

I am trying to fit the parameters to a certain dataset.

LsqFit does not seem to like multivalues function but maybe I am not using it correctly:

```
using CSV
data = CSV.read("panel.csv")
N = size(data)[1]
function fitmodel(x,p)
r =@view x[:,1]
T1 =@view x[:,2]
θ =@view x[:,3]
α1,α2,k,T2 = p
vcat( α1.*T1.^4 .+ k.*(T1.-T2) .- cosd.(θ)./r.^2 , α2.*T2^4 .+ k .* (T2.-T1) )
end
x = convert(Matrix,data[:,[:Distance,:Temp,:Tilt]])
y = zeros(2*N)
# initial guess
r, T1, θ = x[end,:]
α1 = (cosd(θ)/r^2)/T1^4
α2 = α1
dT = 10.0
T2 = T1 - dT
k = α2*T2^4 / dT
p0 = [ α1 α2 k T2 ]
println("Initial guess: ",p0)
using LsqFit
fit = curve_fit(fitmodel, x, y, p0)
```

I get the following error:

```
ERROR: MethodError: no method matching levenberg_marquardt(::NLSolversBase.OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,2}}, ::Array{Float64,2})
You might have used a 2d row vector where a 1d column vector was required.
Note the difference between 1d column vector [1,2,3] and 2d row vector [1 2 3].
You can convert to a column vector with the vec() function.
Closest candidates are:
levenberg_marquardt(::NLSolversBase.OnceDifferentiable, ::AbstractArray{T,1}; x_tol, g_tol, maxIter, lambda, tau, lambda_increase, lambda_decrease, min_step_quality, good_step_quality, show_trace, lower, upper, avv!) where T at /home/bgodard/.julia/packages/LsqFit/LeVh7/src/levenberg_marquardt.jl:34
Stacktrace:
[1] lmfit(::NLSolversBase.OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,2}}, ::Array{Float64,2}, ::Array{Float64,1}; autodiff::Symbol, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/bgodard/.julia/packages/LsqFit/LeVh7/src/curve_fit.jl:68
[2] lmfit(::NLSolversBase.OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,2}}, ::Array{Float64,2}, ::Array{Float64,1}) at /home/bgodard/.julia/packages/LsqFit/LeVh7/src/curve_fit.jl:68
[3] lmfit(::LsqFit.var"#18#20"{typeof(fitmodel),Array{Float64,2},Array{Float64,1}}, ::Array{Float64,2}, ::Array{Float64,1}; autodiff::Symbol, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/bgodard/.julia/packages/LsqFit/LeVh7/src/curve_fit.jl:64
[4] lmfit(::Function, ::Array{Float64,2}, ::Array{Float64,1}) at /home/bgodard/.julia/packages/LsqFit/LeVh7/src/curve_fit.jl:61
[5] curve_fit(::typeof(fitmodel), ::Array{Float64,2}, ::Array{Float64,1}, ::Array{Float64,2}; inplace::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/bgodard/.julia/packages/LsqFit/LeVh7/src/curve_fit.jl:115
[6] curve_fit(::Function, ::Array{Float64,2}, ::Array{Float64,1}, ::Array{Float64,2}) at /home/bgodard/.julia/packages/LsqFit/LeVh7/src/curve_fit.jl:106
```

Both the model and observation are column vectors but obviously I need to have the observation and modeled vectors twice as long as the input data.

I have also tried an observation matrix with:

```
y = zeros(N,2)
```

and using hcat instead of vcat in the model return value, but this results in the same error.

Here I am targeting (0,0) for the observation so I could create a univalue function returning the vector norm but this would not be numerically optimal. Another solution would be to duplicate the data vector entries and add a discrete independent variable to switch between the 2 components of the function but this looks convoluted.

What is the easiest way with the Julia ecosystem to fit a multi-value function?