For interpolation linear solves, don't loop over npt Lagrange polynomials - scipy.linalg.solve_triangular can handle multiple RHS (https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve_triangular.html) Should lead to a speed improvement.