In this author's opinion, the most exciting recent prospect is the
application of rescaling of time discussed above, particularly by
Coomes et al. [18]. The fact that they could
shadow the Lorenz equations 
times longer using a rescaling
of time is remarkable. Miller [58] has a diagram showing
the distance between two nearby solutions with 8 particles, in which
2 of the particles form a temporary binary star with a highly
elliptical orbit. There are very large spikes in the separation of
the two solutions,
representing the times of closest approach of the two
stars, because at closest approach the stars speed up immensely.
Thus, the phase space separation of the two nearby solutions temporarily
increases dramatically. Without a re-scaling of time, it is unlikely
that a shadowing algorithm would be able to shadow such an orbit,
when in fact it is not unreasonable to say that a true solution that
temporarily moves away from a numerical solution in this fashion can
still be labelled a shadow. It is already known that glitches in
N-body simulations occur most frequently near close approaches
[62, 63, 33]. It will be interesting to see how often
the glitches can be eliminated by a rescaling of time. It may be
related to Chow & Palmer's 
-- perhaps p measures the
number of timesteps in an encounter.
Another interesting point in regard to non-hyperbolicity is the
observation that N-body systems have some purely imaginary eigenvalues
(cf. footnote here 
). This is not surprising,
since there are components of regular motion in most N-body systems.
It may be that these ``rotating''
directions can not be ``pinned to zero'' at the endpoints of the
shadow. Perhaps it may be possible to handle this non-hyperbolic
structure by allowing these ``rotating'' directions to have non-zero
boundary conditions at the endpoints, unlike the exponential directions.
Quinlan and Tremaine [62, 63] ask the question of whether shadows, even if they exist, are typical of true orbits chosen at random. It would seem odd if they were not. Perhaps bi-shadowing may be applicable. Bi-shadowing studies not only the existence of a true solution close to a numerical solution, but the existence of possible numerical solutions close to arbitrary true solutions.
Many authors have demonstrated that the rate of exponential divergence decreases as a system becomes less collisional [47, 48, 49, 27, 63]. Hayes [33] demonstrates that a simple scaling of the N-body problem gives a shadow length inversely proportional to N, which seems rather disappointing. However, Quinlan and Tremaine, and Hayes, discuss the possibility that, even if shadowing becomes more difficult in higher dimensions, perhaps the fact that N-body systems become smoother and less collisional with higher N may compensate for this, allowing smooth systems with large numbers of particles to be shadowable.
There are many, many more interesting areas of further research. The author's Master's Thesis [33] discusses many of them. It also discusses other measures of error that may be used if shadowing turns out to be too stringent a measure of error.
Finally, there may be the possibility of proving rigorously that a shadow does not exist. The literature on this subject seems rather sparse; it is mentioned briefly in Adams et al. [2] and Iserles et al. [42].