NDSolve and a vector function [duplicate]

Question

I'm trying to solve a set of differential equations which include a numerical function that takes a vector and returns a vector.

f[x_?NumericQ, y_?NumericQ] := {x^2 + y^2, x - y}

When I try and use this in NDSolve, Mathematica can't tell that it is a vector and returns NDSolve::underdet.

NDSolve[{x[0] == 0, y[0] == 1, x'[t] + Sin[y[t]] + z1[t] == 0, 
  y'[t] == -Cos[x[t]] + z2[t], {z1[t], z2[t]} == f[x[t], y[t]]}, {x,y, z1, z2}, {t, 0, 1}]

I can work around this by splitting f into parts, as follows, but that seems like it may be bad for performance (and looks really inelegant).

fx[x_?NumericQ, y_?NumericQ] := f[x, y][[1]]
fy[x_?NumericQ, y_?NumericQ] := f[x, y][[2]]

So this works:

NDSolve[{x[0] == 0, y[0] == 1, x'[t] + Sin[y[t]] + z1[t] == 0, 
  y'[t] == -Cos[x[t]] + z2[t], z1[t] == fx[x[t], y[t]], z2[t] == fy[x[t], y[t]]}, 
  {x, y, z1, z2}, {t, 0, 1}]

but that requires me to memo-ize f and define some kind of i-th part that looks really clunky when f is bigger. My actual application has $n$ different equations for $4*n$ unknowns, plus $n$ of these equations representing forcing terms, where $n$ is a number of points of spatial discretization of a PDE, which is currently 30.

Answer 1

You could just use Thread to turn an equality between 2 vectors in to a list of equalities,

Thread[{r, b, c} == {1, 2, 3}]
(* {r == 1, b == 2, c == 3} *)

Of course, in your case, since you have NumericQ test on f, this doesn't work out of the box. You can get around that with With,

f[x_, y_] := {x^2 + y^2, x - y};
With[
 {f2 = f[x[t], y[t]]},
 NDSolve[{x[0] == 0, y[0] == 1, x'[t] + Sin[y[t]] + z1[t] == 0, 
   y'[t] == -Cos[x[t]] + z2[t], Thread[{z1[t], z2[t]} == f2]}, {x, y, 
   z1, z2}, {t, 0, 1}]
 ]