Make solve() able to return newton stages

[niklasreich]

Currently it is not possible for the solve() command to return the stages of the Newton scheme. This information is needed for the EIM in nonlinear problems, so for now one has to use the newton() command directly. It would be nice to be able to make the solve() command return the Newton stages as well.

[sdrave]

This is a bit more difficult than I initially thought: While returning the Newton stages via the compute output dict is no problem, StationaryModel._compute_solution calls operator.apply_inverse, which then dispatches into newton. Thus, to get hold of the stages, we would need to extend the apply_inverse interface to also let it return additional data. I'm not sure at this point if we want to go this way.