Indexed and array-valued interpolations

Indexed is hardly well-documented, so it is hard to say what is or what is not a bug. However, I will argue that what you encountered is not a bug. The basis of my argument is that what documentation there is indicates the first argument given to Indexed should be an expression that is a valid first argument for Part. An interpolation function does not meet this condition.

f = NDSolve[{y'[t] == y[t], y[0] == {1, 2}}, y, {t, 0, 1}][[1, 1, 2]]
Table[Part[f[t], i], {i, 2}]

Part::partw: Part 2 of InterpolatingFunction[{{0.,1.}},{5,3,1,{26},{4},0,0,0,0,Automatic,{},{},False},{{<<1>>}},{{{1.,2.},{1.,2.}},<<24>>,{{2.71828,5.43656},{2.71828,<<19>>}}},{Automatic}][t] does not exist. >>

I am arguing that one can only use Indexed safely in situations where

 Indexed`[<expression>, <index>] /. Indexed -> Part

produces a valid result.

[Ref: Please see The clearest way to represent Mathematica's evaluation sequence for WReach's traceView2 function, which I will use to trace evaluation.]

I'm a little more confused than I like to be when I answer. Perhaps sharing what I think I know will prompt others who better understand evaluation in Mathematica to contribute.

One of the mysterious things about Mathematica is how it determines when to evaluate an expression as it moves through the digestion process of evaluation. No doubt it's perfectly clear to those who know the internal workings. What I mean is that, putting held expressions and so forth aside, once Mathematica has evaluated an InterpolatingFunction or an Integrate that results in the same (or even a different) expression, it does not need to reevaluate it -- unless something changes which either does or might lead to a different result. There must be something like a been-there-done-that flag associated with the expression, but however it is done, I have to hope and treat it like it does not matter, because I don't how it's done.

Working answer: What seems to be happening is that inside Indexed something happens that suggests to Mathematica that something changed that might lead to a different result for the InterpolatingFunction. The InterpolatingFunction is reevaluated and that causes Indexed to be reevaluated. This leads to an infinite loop, since inside Indexed, once again the InterpolatingFunction is reevaluated, etc.

If you wish to explore this, here is a minimal example:

if = Interpolation[Range@2, InterpolationOrder -> 1];
Indexed[if, 2]

Unlike the example in the OP, using Indexed this way does not make sense; however, it does give the same $IterationLimit::itlim behavior as the original with a minimal-sized trace.

You might first want to see the evaluation process that constructs an InterpolatingFunction. First Interpolation constructs an InterpolatingFunction expression: The standard evaluation sequence, head Interpolation, first argument Range@2, etc. results in

InterpolatingFunction[{{1, 2}}, {5, 3, 0, {2}, {2}, 0, 0, 0, 0, 
  Automatic, {}, {}, False}, {{1, 2}}, {{1}, {2}}, {Automatic}]

which is itself evaluated in six steps.

traceView2[Interpolation[Range@2, InterpolationOrder -> 1]]

The final result has the same input form, but is formatted(!). I don't know for sure why it is formatted, but it seems consistent that it is formatted when Mathematica is done evaluating the expression. I have assumed this and used it as a key in interpreting the trace of the evaluation.

The evaluation of InterpolatingFunction[..] also follows the standard evaluation sequence and results in an expression with the same input form, although now that it is evaluated, it will be formatted when returned as seen above.

Another trick to help read the trace output is to replace all InterpolationFunction[..] expressions by "IFN". You lose the formatting, but the output is smaller and easier to deal with. Since the expressions in a trace are wrapped in HoldForm, this replacement does not cause reevaluation of the trace output.

Here we can see the start of the infinite loop described in the working answer. After if is evaluated, Indexed[IFN, 2] evaluates to Indexed[IFN, {2}], which begins the loop. Note that the InterpolatingFunction is reevaluated (the 7-step block, not expanded). So after the head Indexed and the arguments IFN and {2}, Indexed[IFN, {2}] is evaluated. Instead of evaluation stopping at this point, as it should, the evaluation process is repeated.

Block[{$IterationLimit = 20},
  traceView2[Indexed[if, 2]]
  ] /. i_InterpolatingFunction :> "IFN"

$IterationLimit::itlim: Iteration limit of 20 exceeded. >>

I haven't found another function with this problem, so far. Only 5000 more to check. :)

This is an extended comment, not an answer.

After running into problems similar to yours with evaluations that returned rules involving interpolation functions, I decided to cut the Gordian knot and simply extract the interpolation function(s) for use in any further computations. I have never regretted adopting this strategy and it certainly works well for this problem.

f = NDSolve[{y'[t] == y[t], y[0] == {1, 2}}, y, {t, 0, 1}][[1, 1, 2]];
Plot[f[t], {t, 0., 1.}]

Voila, no fuss, no muss, no bother (old tv commercial catch phrase).

Update

Addressing an issue raised in a comment.

{f1, f2} = Table[With[{i = i}, f[#][[i]] &], {i, 2}];
Plot[{f1[t], f2[t]}, {t, 0., 1.}]