Wrapper for x2sys_solve

**Description of the desired feature**

Implement [`x2sys_solve`](https://docs.generic-mapping-tools.org/6.2/supplements/x2sys/x2sys_solve.html) which "Determine least-squares systematic correction from crossovers". Part of the ['x2sys`](https://docs.generic-mapping-tools.org/6.2/modules.html#x2sys) suite of modules!

This is a feature requested at https://forum.generic-mapping-tools.org/t/help-with-pygmt-x2sys-cross/1931/6 by [Thor524](https://forum.generic-mapping-tools.org/u/thor524), and will be a follow up of #545. Note that [`x2sys_list`](https://docs.generic-mapping-tools.org/6.2/supplements/x2sys/x2sys_list.html) "Extract subset from crossover data base" may also need to be implemented for `x2sys_solve` to work. The workflow is something as follows:

`x2sys_init` :heavy_check_mark: -> `x2sys_cross` :heavy_check_mark: -> `x2sys_list` :negative_squared_cross_mark: -> `x2sys_solve` :negative_squared_cross_mark:

Reference:

- Wessel, P. (2010). Tools for analyzing intersecting tracks: The x2sys package. Computers & Geosciences, 36(3), 348–354. https://doi.org/10.1016/j.cageo.2009.05.009

**Are you willing to help implement and maintain this feature?** Not on the top of my priority, but anyone is free to give this a try, following the previous `x2sys_init/cross` Pull Request at #546. Let me or the PyGMT team know if you're interested in helping and we can guide you along.