code: def maxflow(V,M,source,sink):
      
      # create max-flow underlying model, on which to find cuts
      model = Model("maxflow")
      
      f = {} # flow variable
      for (i,j) in M:
          f[i,j] = model.addVar(lb=-M[i,j], ub=M[i,j], name="flow(%s,%s)"%(i,j))
  
      cons = {}
      for i in V:
          if i != source and i != sink:
              cons[i] = model.addCons(
                  quicksum(f[i,j] for j in V if i<j and (i,j) in M) - \
                  quicksum(f[j,i] for j in V if i>j and (j,i) in M) == 0,
                  "FlowCons(%s)"%i)
  
      model.setObjective(quicksum(f[i,j] for (i,j) in M if i==source), "maximize")
  
      # model.write("tmp.lp")
      model.data = f,cons
      return model
docs: maxflow: maximize flow from source to sink, taking into account arc capacities M
      Parameters:
          - V: set of vertices
          - M[i,j]: dictionary or capacity for arcs (i,j)
          - source: flow origin
          - sink: flow target
      Returns a model, ready to be solved.