Match the precision level in the lens equation solver between time delays and image positions

[jiwoncpark]

When combining 100s of simulated PEMD lenses, the H0 recovery precision (0.4%) approaches a point where the numerical precision level of the lens equation solver starts to matter. My solution for forcing the precision level to be the same between the image positions and time delays has simply been to have the TDCosmology's time_delays method return the image positions that were used to evaluate the time delays. But this will most likely break people's code -- this issue is meant to be a discussion more than a simple code fix!

[sibirrer]

@jiwoncpark I updated the design such that all the lens equation solver specific arguments are part of kwargs_lens_eqn_solver within the kwargs_model arguments and those should be able to be passed to all relevant classes self-consistently, including the ones in the Analysis module. This should also make PR #185 redundant. Let me know if you have questions and whether you want to keep this issue open. Thanks