code: def d2Ibr_dV2(Ybr, V, lam):
      
      nb = len(V)
      diaginvVm = spdiag(div(matrix(1.0, (nb, 1)), abs(V)))
  
      Haa = spdiag(mul(-(Ybr.T * lam), V))
      Hva = -1j * Haa * diaginvVm
      Hav = Hva
      Hvv = spmatrix([], [], [], (nb, nb))
  
      return Haa, Hav, Hva, Hvv
docs: Computes 2nd derivatives of complex branch current w.r.t. voltage.