How about adding `solve!` and also a default implementation for`solve`? ```julia solve(args...; kwargs...) = solve!(init(args...; kwargs...)) ```