Vector-valued ParametricNDSolve, solving for a combination

Question

I am integrating two vector differential equations using ParametricNDSolve, one for $\mathbf{Y}$ and one for $\mathbf{Z}$, and then I'm interested in a combination of the two of them. ParametricNDSolve allows for specifying the output to just be a function of the two, for instance $\mathbf{Y} \cdot \mathbf{Z}$, but the combination I want is (for the case of sixth order vector equations): $$\psi = y_1 z_6 -y_2 z_5 + y_3 z_4 + y_4 z_3 - y_5 z_2 +y_6 z_1$$

Is there a way to get ParametricNDSolve to give me this combination directly?

Example code:

A = {{0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}, {q, 0, -Sin[x], 0}};
Q[x_, q_] = {{0, 1, 0, 0, 0, 0}, {0, 0, 1, 1, 0, 0}, {0, -Sin[x],0, 0, 1, 0}, 
           {0, 0, 0, 0, 1, 0}, {-q, 0, 0, -Sin[x], 0, 1}, {0, -q, 0, 0, 0, 0}};
n = Length[Q[x, q]];
cA = With[{code = N@A}, Compile[{{x, _Real}, {q, _Real}}, code]];
F[x_?NumericQ, q_?NumericQ, Y_?VectorQ] := 
      (Q[x, q] - DiagonalMatrix[ConstantArray[1, n]]).Y;
xa = -4; xb = 4; xm = 0; 
k = -((-xb + xm)/(xa - xm)); c = -((-xa xm + xb xm)/(xa - xm));
Y0 = {0, 0, 0, 0, -1, 0};
Z0 = {0, 1, 0, 0, 0, 0};

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];

I'm interested in finding where $\psi=0$, which I can do with this construction, but is there a more direct way to get it out of the ParametricNDSolveValue?

psi[q_?NumericQ] := (Y[1] Z[6] - Y[2] Z[5] + Y[3] Z[4] + Y[4] Z[3] - 
     Y[5] Z[2] + Y[6] Z[1]) /. {Y[a_] :> YZsol[q][[1]][xm][[a]], 
    Z[a_] :> YZsol[q][[2]][xm][[a]]};

FindRoot[psi[q], {q, 3}]

Edit:

Carl's answer does exactly what I need, but I see a strange timing effect, where after I use FindRoot it takes 6 times longer to do the same calculation:

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]] 
(* {0.259752, Null} *)

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]]
(* {1.69719, Null} *)

Answer 1

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).

Answer 2

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}} *)