# Question about using NDSolve with a matrix equation

I tried to use NDSolve to study a matrix-valued differential backward equation. The following program seems to work:

    n = 15;
tin = 0;
tfin  = 1;

A[t_] := t/5;
Abarra[t_] := -t/2;
Enne[t_] := t + 1;
q[t_] := 2 t/10 tfin;
qbarra[t_] := 3 t/10 tfin;
sbarra[t_] := t/10 tfin + 1;
G[t_] := t/2 ;
M[t_] := t/7;
(*Definition of the involved matrices*)
Z[t_] := 1/(n - 1) Table[
If[i == j && i == 1, A[t] (n - 1),
If[i == j && i > 1, (A[t] + G[t]) (n - 1),
If[j > i && i == 1, Abarra[t], Abarra[t] + M[t]]
]
]
, {i, n}, {j, n}];
Q[t_] := Table[
If[i == j == 1, q[t] + qbarra[t],
If[i == 1 || j == 1, -1/(n - 1) sbarra[t] qbarra[t],
1/(n - 1)^2 sbarra[t] qbarra[t] sbarra[t]
]
]
, {i, n}, {j, n}]
Mim[t_] := Table[If[i == j == 1, Enne[t]^(-1), 0], {i, n}, {j, n}];
Zero = Table[0, {i, n}, {j, n}];

func[X_?MatrixQ, t_] :=  Z[t].X + X.Transpose[Z[t]] - X.Mim[t].X + Q[t];
sol = NDSolve[{P'[t] == - func[P[t], t], P[tfin] == Q[tfin]}, P, {t, tfin, tin}]


I have now the following questions:

1. How do I plot the graphic of P's entries? I tried to write

Plot[Evaluate[P[t][[1, 2]] /. sol], {t, 0, tfin}]


but that gives an error and moreover is different from what I expect because of the final conditions.

2. If I change a bit the NDSolve argument as

sol = NDSolve[{P'[t] + func[P[t], t] == Zero, P[tfin] == Q[tfin]}, P, {t, tfin, tin}]


then I have an error. Why?

For 1, the problem is that you want to evaluate the ReplaceAll (/.) but not evaluate Part ([[1, 2]]) until a matrix value of the solution is computed. The evaluation has to be separated. One way is to use With:

With[{pfn = P /. First[sol]},
Plot[pfn[t][[1, 2]], {t, 0, tfin}]
] I cannot say whether it is what you expect, obviously. :)

For 2, I haven't fully tracked down the problem. I don't get an error in V10, but I get one in V9 that suggests that it could have to do with threading. Edit: This now seems likely after examining the Trace of the command in V9, which is reflected in the following. Evaluating the solution for P'[t], that is,

P'[t] -> - func[P[t], t] + Zero


results in a 15 x 15 matrix, in which each entry is -func[P[t], t]. In V10, NDSolve could be solving 15^2 DE's, each of which is a 15 x 15 system. It runs for a while -- I haven't let it go more than a few minutes, though. In V9, NDSolve might miscompute the dimension, which is the error it gives:

NDSolve::ndfdmc: Computed derivatives do not have dimensionality consistent with the initial conditions. >>

Aside from not using Zero, another workaround is to define

Zero[t_?NumericQ] = ConstantArray[0, {n, n}];


And then change Zero to Zero[t]:

sol = NDSolve[{P'[t] + func[P[t], t] == Zero[t], P[tfin] == Q[tfin]}, P, {t, tfin, tin}]


Then Zero[t] will not evaluate to a matrix until the other elements in the equation evaluate to matrices.

P.S. Generally, it is advisable to avoid starting your own variables and functions with capital letters. That way they will never conflict with built-in variables and functions.