When given the task to change my explicit ocean model to an implicit version, I had assumed that I would need to invert the ‘A’ matrix to solve Ax=B, and then do some matrix math. I found a few web sites that warn against that because it is possibly a lot more computationally intensive than solving directly. John Cook’s blog explains more about why.

Python’s NumPy has linalg.solve(A, B), which returns the ‘x’ array

x = numpy.linalg.solve(A,B)

It uses a LU decomposition method for solving (not inversion).