Why I can not get the plot when I use NDSolve`ProcessSolutions?

Question

Why I can not get the plot when I use NDSolve`ProcessSolutions? Anyone can give me a clue. Thanks a lot!

state = First[NDSolve`ProcessEquations[{
               D[u[t, x], t] == D[u[t, x], x, x],
               u[0, x] == Cos[π /2 x], u[t, 0] == 1, u[t, 1] == 0},
               u, t, {x, 0, 1}]];
NDSolve`Iterate[state, {0, 4}];
NDSolve`ProcessSolutions[state, "Forward"];
Plot[Evaluate[u[4, x] /. %], {x, 0, 1}]

Answer 1

Actually a very simple modification would make your approach work:

Plot[Evaluate[u[4., x] /. %], {x, 0, 1}]

the problem is that for pattern matching (which is what ReplaceAll (/.) does in the background) 4 (with head Integer) and 4. (with head Real) are two different things. You even might meet cases where using 4. would not be enough to make this work, e.g. when the end point of the integration for some reason is not exactly four (I'm not sure whether or not that could happen with NDSolve). So in general pattern matching using (machine precision) reals is error prone. There is actually a very simple approach to extract the solution functions, but it might fail in similar cases, you have to adopt the indexing according to the expression from which you are extracting:

state = First[
    NDSolve`ProcessEquations[{
     D[u[t, x], t] == D[u[t, x], x, x],
     u[0, x] == Cos[Pi/2 x],
     u[t, 0] == 1,
     u[t, 1] == 0
    }, u, t, {x, 0, 1}]
];
NDSolve`Iterate[state, {0, 4}];
solutions = NDSolve`ProcessSolutions[state, "Forward"];
usol = solutions[[1, 2, 0]]
dusol = solutions[[2, 2, 0]]

and you then can easily plot the results with e.g.:

Plot[usol[x], {x, 0, 1}]

If you want a solution which will not rely on the exact order and structure of the result (similar to what your original pattern matching solution did) but want the pattern matching to not rely on real numbers to be SameQ, you can reverse the role of expression (u[4.,x]) and pattern (the LHS of the rules in solutions). That can be done using Cases with an appropriate pattern which unfortunately becomes relatively complicated for this case:

 usol = First[Cases[
    solutions, Verbatim[Rule][u[_, _], (i_InterpolatingFunction)[_]]] :> i
 ]]

This makes use of the fact that with Cases you can not only extract subexpressions which match a pattern but also replace those in one go if the second argument is a Rule or RuleDelayed. Here is the same thing if you want to extract the derivative:

dusol = First[Cases[solutions,
    Verbatim[Rule][Derivative[1, 0][u][_, _],(i_InterpolatingFunction)[_]] :> i
]]

Here is another version which makes use of a somewhat simpler pattern and extracts all solutions (but it relies on the fact that ProcessSolutions does return its result in a certain order):

{usol, dusol} = Cases[solutions, (i_InterpolatingFunction)[x] :> i, {2}]

The third argument for Cases is the level(s) in which you want it to search for subexpressions that match the given pattern, in this case that can to be set to only level 2 as that is the level where the interpolating functions sit, but you could also use Infinity which would search in all levels (that last case was my initial answer). You can find all these admittedly nontrivial details to be documented in the online documentation e.g. PatternsAndTransformationRules and Cases.