code: def solve(preconditioning=True):
      """Compute the solution"""
      count = [0]
  
      def residual(P):
          count[0] += 1
  
          d2x = zeros_like(P)
          d2y = zeros_like(P)
  
          d2x[1:-1] = (P[2:] - 2 * P[1:-1] + P[:-2]) / hx / hx
          d2x[0] = (P[1] - 2 * P[0] + P_left) / hx / hx
          d2x[-1] = (P_right - 2 * P[-1] + P[-2]) / hx / hx
  
          d2y[:, 1:-1] = (
              (P[:, 2:] - 2 * P[:, 1:-1] + P[:, :-2]) / hy / hy
          )
          d2y[:, 0] = (P[:, 1] - 2 * P[:, 0] + P_bottom) / hy / hy
          d2y[:, -1] = (P_top - 2 * P[:, -1] + P[:, -2]) / hy / hy
  
          return d2x + d2y + 5 * cosh(P).mean() ** 2
  
      # preconditioner
      if preconditioning:
          M = get_preconditioner()
      else:
          M = None
  
      # solve
      guess = zeros((nx, ny), float)
  
      sol = root(
          residual,
          guess,
          method="krylov",
          options={"disp": True, "jac_options": {"inner_M": M}},
      )
      print("Residual", abs(residual(sol.x)).max())
      print("Evaluations", count[0])
  
      return sol.x
  
imports: from scipy.optimize import root
  from scipy.sparse import spdiags, kron
  from scipy.sparse.linalg import spilu, LinearOperator
  from numpy import cosh, zeros_like, mgrid, zeros, eye
label_desc: Use image node_id 2 for calling a global function with example usage: solve(preconditioning) and returns: sol, unknown