docs: def stringfy(expr, sym_const=None, sym_states=None, sym_algebs=None):
      """Convert the right-hand-side of an equation into CVXOPT matrix operations"""
      if not sym_const:
          sym_const = []
      if not sym_states:
          sym_states = []
      if not sym_algebs:
          sym_algebs = []
      expr_str = []
      if type(expr) in (int, float):
          return expr
      if expr.is_Atom:
          if expr in sym_const:
              expr_str = 'self.{}'.format(expr)
          elif expr in sym_states:
              expr_str = 'dae.x[self.{}]'.format(expr)
          elif expr in sym_algebs:
              expr_str = 'dae.y[self.{}]'.format(expr)
          elif expr.is_Number:
              if expr.is_Integer:
                  expr_str = str(int(expr))
              else:
                  expr_str = str(float(expr))
              # if expr.is_negative:
              #     expr_str = '{}'.format(expr)
              # else:
              #     expr_str = str(expr)
          else:
              raise AttributeError('Unknown free symbol <{}>'.format(expr))
      else:
          nargs = len(expr.args)
          arg_str = []
          for arg in expr.args:
              arg_str.append(stringfy(arg, sym_const, sym_states, sym_algebs))
  
          if expr.is_Add:
              expr_str = ''
              for idx, item in enumerate(arg_str):
                  if idx == 0:
                      if len(item) > 1 and item[1] == ' ':
                          item = item[0] + item[2:]
                  if idx > 0:
                      if item[0] == '-':
                          item = ' ' + item
                      else:
                          item = ' + ' + item
                  expr_str += item
  
          elif expr.is_Mul:
              if nargs == 2 and expr.args[0].is_Integer:  # number * matrix
                  if expr.args[0].is_positive:
                      expr_str = '{}*{}'.format(*arg_str)
                  elif expr.args[0] == Integer('-1'):
                      expr_str = '- {}'.format(arg_str[1])
                  else:  # negative but not -1
                      expr_str = '{}*{}'.format(*arg_str)
              else:  # matrix dot multiplication
                  if expr.args[0] == Integer('-1'):
                      # bring '-' out of mul()
                      expr_str = ', '.join(arg_str[1:])
                      expr_str = '- mul(' + expr_str + ')'
                  else:
                      expr_str = ', '.join(arg_str)
                      expr_str = 'mul(' + expr_str + ')'
          elif expr.is_Function:
              expr_str = ', '.join(arg_str)
              expr_str = str(expr.func) + '(' + expr_str + ')'
          elif expr.is_Pow:
              if arg_str[1] == '-1':
                  expr_str = 'div(1, {})'.format(arg_str[0])
              else:
                  expr_str = '({})**{}'.format(*arg_str)
          elif expr.is_Div:
              expr_str = ', '.join(arg_str)
              expr_str = 'div(' + expr_str + ')'
          else:
              raise NotImplementedError
      return expr_str