Vector-valued ParametricNDSolve, solving for a combination

Define a quadratic form matrix:

form = {
    {0, 0, 0, 0, 0, 1},
    {0, 0, 0, 0, -1, 0},
    {0, 0, 0, 1, 0, 0},
    {0, 0, 1, 0, 0, 0},
    {0, -1, 0, 0, 0, 0},
    {1, 0, 0, 0, 0, 0}
};

Then:

YZsol2 = ParametricNDSolveValue[
    {Y'[x] == F[x,q,Y[x]], Y[xa]==Y0,
    Z'[x] == 1/k F[k x+c,q,Z[x]], Z[xa]==Z0},
    Y[xm] . form . Z[xm],
    {x,xa,xm},
    q
];

will produce the desired function. For example:

YZsol2 /@ Range[0, 4]
psi /@ Range[0, 4]

{-24.1853, -148.784, -6240.84, -7841.06, 12815.}

{-24.1853, -148.784, -6240.84, -7841.06, 12815.}

Another possibility is to include your desired output as an additional equation (constraint):

YZsol3 = ParametricNDSolveValue[
    {
    Y'[x] == F[x,q,Y[x]], Y[xa]==Y0,
    Z'[x] == 1/k F[k x+c,q,Z[x]], Z[xa]==Z0,
    U[x] == Y[x] . form . Z[x]
    },
    U[xm],
    {x,xa,xm},
    q
];

Check:

YZsol3 /@ Range[0, 4]

{-24.1853, -148.784, -6240.84, -7841.06, 12815.}

Note that there is a bug as of June 2019 with using this inside a FindRoot, which makes it much slower unless we add Method -> {"ParametricSensitivity" -> None} as detailed in Vector ParametricNDSolve and FindRoot interaction).

A hybrid of yours and @CarlWoll's approach seems to perform well after FindRoot.

Clear[YZsol, psi2]
form = Reverse@DiagonalMatrix@PadRight[{}, 6, {1, -1, 1}];
YZsol = ParametricNDSolveValue[{Y'[x] == F[x, q, Y[x]], Y[xa] == Y0, 
    Z'[x] == 1/k F[k x + c, q, Z[x]], Z[xa] == Z0}, {Y, Z}, {x, xa, 
    xm}, q];
psi2[q_?NumericQ] := (Y.form.Z) /. {Y :> YZsol[q][[1]][xm], 
    Z :> YZsol[q][[2]][xm]};
FindRoot[psi2[q], {q, 3}] 
AbsoluteTiming[psi2 /@ Range[20]]

Clear[YZsol2];
YZsol2 = ParametricNDSolveValue[{Y'[x] == F[x, q, Y[x]], Y[xa] == Y0, 
    Z'[x] == 1/k F[k x + c, q, Z[x]], Z[xa] == Z0}, 
   Y[xm].form.Z[xm], {x, xa, xm}, q];
AbsoluteTiming[YZsol2 /@ Range[20]]


Clear[YZsol3]; YZsol3 = 
 ParametricNDSolveValue[{Y'[x] == F[x, q, Y[x]], Y[xa] == Y0, 
   Z'[x] == 1/k F[k x + c, q, Z[x]], Z[xa] == Z0}, 
  Y[xm].form.Z[xm], {x, xa, xm}, q];
FindRoot[YZsol3[q], {q, 3}];
AbsoluteTiming[YZsol3 /@ Range[20]]
(* {q -> 3.55393} *)
(* {0.279002, {-148.784, -6240.84, -7841.06, 12815., 63393.7, 
  134180., 197976., 214749., 
  139054., -71335.6, -446840., -998422., -1.71155*10^6, -2.5422*10^6, \
-3.41516*10^6, -4.22468*10^6, -4.83742*10^6, -5.09767*10^6, \
-4.83432*10^6, -3.86961*10^6}} *)
(* {0.337883, {-148.784, -6240.84, -7841.06, 12815., 63393.7, 
  134180., 197976., 214749., 
  139054., -71335.6, -446840., -998422., -1.71155*10^6, -2.5422*10^6, \
-3.41516*10^6, -4.22468*10^6, -4.83742*10^6, -5.09767*10^6, \
-4.83432*10^6, -3.86961*10^6}} *)
(* {2.20664, {-148.785, -6240.83, -7841.05, 12815., 63393.6, 
  134180., 197976., 214748., 
  139054., -71335.6, -446840., -998421., -1.71155*10^6, -2.5422*10^6, \
-3.41516*10^6, -4.22467*10^6, -4.83742*10^6, -5.09767*10^6, \
-4.83432*10^6, -3.86961*10^6}} *)