def matrix_free_solve(matvec, b):
    return lax.custom_linear_solve(matvec, b, explicit_jacobian_solve, explicit_jacobian_solve)
def matrix_free_solve(matvec, b):
    return lax.custom_linear_solve(matvec, b, explicit_jacobian_solve, explicit_jacobian_solve)