code: def energy_corrections(perturbation, n, *, a=10, mass=0.5):
      """
      Calculating first two order corrections due to perturbation theory and
      returns tuple where zero element is unperturbated energy, and two second
      is corrections
  
      ``n``
          the "nodal" quantum number. Corresponds to the number of nodes in the
          wavefunction.  n >= 0
      ``a``
          width of the well. a > 0
      ``mass``
          mass.
  
      """
      x, _a = var("x _a")
  
      Vnm = lambda n, m, a: Integral(
          X_n(n, a, x) * X_n(m, a, x) * perturbation.subs({_a: a}),
          (x, 0, a),
      ).n()
  
      # As we know from theory for V0*r/a we will just V(n, n-1) and V(n, n+1)
      #   wouldn't equals zero
  
      return (
          E_n(n, a, mass).evalf(),
          Vnm(n, n, a).evalf(),
          (
              Vnm(n, n - 1, a) ** 2
              / (E_n(n, a, mass) - E_n(n - 1, a, mass))
              + Vnm(n, n + 1, a) ** 2
              / (E_n(n, a, mass) - E_n(n + 1, a, mass))
          ).evalf(),
      )
  
imports: from sympy.core import pi
  from sympy import Integral, var, S
  from sympy.functions import sin, sqrt
label_desc: Use image node_id 3 for calling a global function with example usage: energy_corrections(perturbation, n) and returns: E_n,  Vnm,  unknown