How to fit parameters in a system of coupled differential equations

I just wanted to point out that all the routes above will work if there is only one differential equation:

data = NDSolveValue[{
     x''[t] - k1*(1 - x[t]^2)*x'[t] + k2*x[t] == 0, 
     x[0] == 2, x'[0] == 0} /. {k1 -> 1., k2 -> 1.}, 
   Table[{t, x[t] + RandomReal[{-.3, .3}]}, {t, 0, 10, .2}], {t, 10}];
dataT = data\[Transpose];
ti = dataT[[1, All]];
di = dataT[[2, All]];
pfun = ParametricNDSolveValue[{
    x''[t] - k1*(1 - x[t]^2)*x'[t] + k2*x[t] == 0,
    x[0] == 2, x'[0] == 0},
    x, {t, 0, 10}, {k1, k2}];

FindFit finds the best-fit parameters sucessfully:

fit = FindFit[data, pfun[k1, k2][t], {{k1, 2}, {k2, 0}}, t]
Out[1]= {k1 -> 1.09028, k2 -> 1.02729}

FindMinimum finds them too:

FindMinimum[Total[(di - Map[pfun[k1, k2], ti])^2], {k1, k2}]
Out[2]= {1.41041, {k1 -> 1.09028, k2 -> 1.02729}}

And the Module approach also produced the same result:

objfun[k_] := Module[{},
  fun = NDSolveValue[{
     x''[t] - k[[1]]*(1 - x[t]^2)*x'[t] + k[[2]]*x[t] == 0,
     x[0] == 2, x'[0] == 0},
    x, {t, 0, 10}];
  Total[(di - Map[fun, ti])^2]
  ]
FindMinimum[objfun[{k1, k2}], {{k1, 1.0}, {k2, 0.0}}]
Out[3]= {1.41041, {k1 -> 1.09028, k2 -> 1.02729}}

Define this function:

f[k1_?NumericQ, k2_?NumericQ, k3_?NumericQ] := 
 Sum[Total[(ci[[i, All]] - Map[pfun[k1, k2, k3][[i]], ti])^2], {i, 1, 3}] // Quiet

Then,

fit = NMinimize[f[k1, k2, k3], {k1, k2, k3}];
params = fit // Last
(*{k1 -> 0.000194805, k2 -> 0.0291469, k3 -> 0.109229}*)

Plot it,

Table[Show[
  ListPlot[Transpose[{ti, ci[[i]]}]],
  Plot[(pfun[k1, k2, k3] /.params)[[i]][t], {t, ti[[1]], ti // Last}
   ]], {i, 1, 3}]