docs: def find_phase_transformation(Ne, Nl, r, Lij, verbose=0,
                                return_equations=False, **kwds):
      """This function returns a phase transformation specified as a list of
      lenght Ne whose elements correspond to each theta_i. Each element is a
      list of length Nl which specifies the coefficients multiplying each
      optical frequency omega^l. So for instance [[1,1],[1,0],[0,0]] means
  
      theta_1=omega^1+omega^2
      theta_2=omega^1
      theta_3=0.
      """
      # We first define the needed variables
      _omega_laser = [Symbol('omega_laser'+str(l+1)) for l in range(Nl)]
      _theta = [Symbol('theta'+str(i+1)) for i in range(Ne)]
  
      # We find all the equations that the specified problem has to fulfil.
      eqs = []
      for i in range(Ne):
          for j in range(i+1, Ne):
              if type(r[0]) == list:
                  if (r[0][i][j] != 0) or (r[1][i][j] != 0) or (r[2][i][j] != 0):
                      for l in range(Nl):
                          if l+1 in Lij[i][j]:
                              eqs += [_omega_laser[l] + _theta[j] - _theta[i]]
              else:
                  if (r[0][i, j] != 0) or (r[1][i, j] != 0) or (r[2][i, j] != 0):
                      for l in range(Nl):
                          if l+1 in Lij[i][j]:
                              eqs += [_omega_laser[l] + _theta[j] - _theta[i]]
  
      if return_equations:
          return eqs
  
      sol = solve(eqs, _theta, dict=True)[0]
      for i in range(Ne):
          if _theta[i] not in sol.keys():
              sol.update({_theta[i]: _theta[i]})
  
      sol_simple = {_theta[i]: sol[_theta[i]]-sol[_theta[-1]] for i in range(Ne)}
  
      sol = []
      for i in range(Ne):
          soli = []
          for l in range(Nl):
              soli += [diff(sol_simple[_theta[i]], _omega_laser[l])]
          sol += [soli]
  
      return sol