Least square with multi-value function

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?

I found the problem. The initial parameter vector was a row. It needs to be a column:

p0 = [ α1, α2, k, T2 ]