# Numerically solving a system of ODEs where the functions are vectorized

I am attempting to solve a system of ODEs where a, b, cv1 and cv2 are vectors.

Tu = {{-3 / 4, 1 / 4}, {5, -6}}
tu = {2/4, 1}
pu = {.7, .3}

a[y_] = {Subscript[a, 1][y], Subscript[a, 2][y]};
b[y_] = {Subscript[b, 1][y], Subscript[b, 2][y]};
cv1[y_] = {Subscript[cv1, 1][y], Subscript[cv1, 2][y]}
cv2[y_] = {Subscript[Subscript[cv2, 1][y], 1], Subscript[Subscript[cv2, 2][y], 2]}

(* Ode system *)
a'[y] == a[y].Tu
b'[y] == Tu.b[y]
cv1'[y] == Tu.cv1[y] + Subscript[a, 1][y].tu
cv2'[y] == Tu.cv2[y] + Subscript[a, 2][y].tu*)

S =
NDSolve[
{a'[y] == a[y].Tu,
b'[y] == Tu.b[y],
cv1'[y] == Tu.cv1[y] + Subscript[a, 1][y].tu,
cv2'[y] == Tu.cv2[y] + Subscript[a, 2][y].tu, a == pu,
b == tu, cv1 == cv2 == 10^-5},
{a[y], b[y], cv1[y], cv2[y] }, {y, 0, ∞}]


I am getting an output informing me that there (incorrectly) exist more dependent variables than equations. How should I specify this system?

### Edit

I fixed the initial code (the above is now the updated version); the previous error no longer occurs, though I now receive the following:

NDSolve::ndnco: The number of constraints (6) (initial conditions) is not equal to the total differential order of the system plus the number of discrete variables (8).

EDIT:

Note that using the system as specified in the answer will not retrieve the correct cv1 and cv2 expressions. Instead specify those in the system as:

cv1'[y] ==
Tu.cv1[y] + (MapIndexed[Indexed[a[y], #2] &, Range] // First)*tu,
cv2'[y] ==
Tu.cv2[y] + (MapIndexed[Indexed[a[y], #2] &, Range] // Last)*tu

• There is no initial data. Jan 22 '19 at 0:01
• My bad, I updated it with the correct initial data and am getting new errors as will be listed Jan 22 '19 at 1:04
• It's just unclear to me what system you're trying to solve. Is Subscript[a, 1][y] and Subscript[a, 2][y] scalar? If so, why there exists Subscript[a, 1][y].tu in the system? (Dot (.) is for products of vectors, matrices, and tensors. ) Is cv2 just a vector like a, b and cv1? If so, why is it defined in such a strange way? Can you show us the system in traditional math notation? Jan 25 '19 at 5:26

Tu = {{-3/4, 1/4}, {5, -6}};

tu = {2/4, 1};
pu = {.7, .3};
ym = 1;

(*Ode system*)
eq = {a'[y] == a[y].Tu, b'[y] == Tu.b[y],
cv1'[y] == Tu.cv1[y] + First[a[y]]*tu,
cv2'[y] == Tu.cv2[y] + Last[a[y]]*tu};
ic = {a == pu, b == tu, cv1 == cv2 == {10^-5, 10^-5}};

S = NDSolve[{eq, ic}, {a, b, cv1, cv2}, {y, 0, ym}]

{Plot[a[y] /. S, {y, 0, ym}, AxesLabel -> {"y", "a"}],
Plot[b[y] /. S, {y, 0, ym}, AxesLabel -> {"y", "b"},
PlotRange -> All],
Plot[cv1[y] /. S, {y, 0, ym}, AxesLabel -> {"y", "cv1"}],
Plot[cv2[y] /. S, {y, 0, ym}, AxesLabel -> {"y", "cv2"}]} If we change the code, excluding First and Last, the result will change, although in symbolic form the equations look the same.

(*Ode system*)
eq = {a'[y] == a[y].Tu, b'[y] == Tu.b[y],
cv1'[y] == Tu.cv1[y] + (a[y].{1, 0})*tu,
cv2'[y] == Tu.cv2[y] + (a[y].{0, 1})*tu};
ic = {a == pu, b == tu, cv1 == cv2 == {10^-5, 10^-5}};
S = NDSolve[{eq, ic}, {a, b, cv1, cv2}, {y, 0, ym}];
{Plot[a[y] /. S, {y, 0, ym}, AxesLabel -> {"y", "a"}],
Plot[b[y] /. S, {y, 0, ym}, AxesLabel -> {"y", "b"},
PlotRange -> All],
Plot[cv1[y] /. S, {y, 0, ym}, AxesLabel -> {"y", "cv1"}],
Plot[cv2[y] /. S, {y, 0, ym}, AxesLabel -> {"y", "cv2"}]} I assume that functions First and Last are ignored when solving equations. @Janeiro offers its solution.

• A colleague who had coded this in R had noted that he had different forms for CV1 and CV2 -- is the First[] and Last[] commands here proper in the system of equations wrt a's elements? Jan 24 '19 at 16:57
• Problem solved. Use : cv1'[y] == Tu.cv1[y] + (MapIndexed[Indexed[a[y], #2] &, Range] // First)*tu, cv2'[y] == Tu.cv2[y] + (MapIndexed[Indexed[a[y], #2] &, Range] // Last)*tu (will update question) Jan 24 '19 at 17:14
• The result is the same. What is the point of adding? Jan 24 '19 at 18:33
• @Janeiro You are right, really NDSolve[]  ignores First[] and Last[]. I published the second solution. Jan 24 '19 at 19:20
• NDSolve doesn't ignore First and Last, because they've already disappeared when the eq=… line is executed: First@a[y] evaluates to y immediately. Notice most functions for list manipulation are not limited to List manipulation e.g. try {First@f[x], Most@f[x,y,z]}`. Jan 25 '19 at 5:08