Linear systems are usually solved in Python with the linalg.solve command of the SciPy package. This function just is a Python interface to the Fortran 90 LAPACK routines which, in turn, are based on the highly optimized BLAS library. Both packages are based upon work supported by the National Science Foundation and the Department of Energy. In the example below, we use this command to solve the BVP of the previous section. Notice that this algorithm is extremely faster than our implementation in the previous section. The reason is that the algorithm is smart and exploits the fact that \(A\) is a "sparse" matrix, namely that most of its elements are zero. This clearly counts a lot when the size of the matrix is large. For instance, for \(N=1000\text{,}\) linalg.solve takes just about 0.1 seconds to solve the system, against the 5 minutes of our implementation of the LU decomposition!!
