diff --git a/examples/mvc/README.md b/examples/mvc/README.md
new file mode 100644
index 0000000000..d4febf6611
--- /dev/null
+++ b/examples/mvc/README.md
@@ -0,0 +1,201 @@
+```python
+import networkx as nx
+import csv
+import matplotlib.pyplot as plt
+
+from mvc import qubo_mvc, qubo_mvc_penalty, mvc_feasibility, mvc_easy_fix, mvc_energy
+```
+
+# MVC Definition
+
+Given an undirected graph with a set of nodes V and edges E, A vertex cover of a graph is a set of nodes that includes at
+least one endpoint of every edge of the graph. The Minimum Vertex Cover (MVC) problem is an optimisation problem
+that find smallest vertex cover of a given graph.
+
+# Load graph from csv
+
+
+```python
+def load_csv(filename:str):
+    """ Load graph from csv file
+    """
+    
+    with open(filename, 'r', newline='') as f:
+        reader = csv.reader(f)
+        data = [[int(row[0]), int(row[1]), float(row[2])] for row in reader]
+    
+    nodes = [k0 for k0, k1, v in data if k0 == k1]
+    edges = [[k0, k1] for k0, k1, v in data if k0 != k1]
+    
+    g = nx.Graph()
+    g.add_nodes_from(nodes)
+    g.add_edges_from(edges)
+    
+    node_weights = {k0: {'weight': v} for k0, k1, v in data if k0 == k1}
+    edge_weights = {(k0, k1): {'weight': v} for k0, k1, v in data if k0 != k1}
+
+    nx.set_edge_attributes(g, edge_weights)
+    nx.set_node_attributes(g, node_weights)
+
+    return g
+```
+
+csv Format:  
+
+0,0,0.4696822179283675  
+1,1,0.6917392308135916  
+2,2,0.7130420958314817  
+3,3,0.49342677332980056  
+4,4,0.49661600571433673  
+5,5,0.55601135361347  
+6,6,0.2831095711244086  
+7,7,0.18737823232636885  
+0,1,0.98602725407379  
+1,2,0.935853268681996  
+2,3,0.23080434023289664  
+4,5,0.25521731195623476  
+5,6,0.37096743935296606  
+6,7,0.04408120693490547  
+0,2,0.5619060764177991  
+4,7,0.48402095290884917  
+1,6,0.7106584134580318  
+
+
+
+```python
+g = load_csv('mvc.csv')
+```
+
+
+```python
+pos = nx.spring_layout(g, seed=1234)
+nx.draw(g, pos)
+```
+
+
+    
+![png](README_files/README_7_0.png)
+    
+
+
+
+```python
+print(f'Number of node is {len(g.nodes)}')
+print(f'Number of edge is {len(g.edges)}')
+```
+
+    Number of node is 8
+    Number of edge is 9
+
+
+# QUBO formulation
+
+## Estimate penalty
+
+
+```python
+penalty = qubo_mvc_penalty(g)
+print(f'The penalty is {penalty}')
+```
+
+    The penalty is 1.6548482220717595
+
+
+## Formulate QUBO (weighted minimum vertex cover)
+
+
+```python
+linear, quadratic = qubo_mvc(g, penalty=penalty)
+print(f'Number of linear terms: {len(linear)}')
+print(f'Number of quadratic terms: {len(quadratic)}\n')
+```
+
+    Number of linear terms: 8
+    Number of quadratic terms: 9
+    
+
+
+# Generate random solutions
+
+
+```python
+import numpy as np
+rng = np.random.default_rng(seed=1234)
+random_solution = {i: rng.integers(2) for i in g.nodes}
+print(f'The random solution is {random_solution}\n')
+```
+
+    The random solution is {0: 1, 1: 1, 2: 1, 3: 0, 4: 0, 5: 1, 6: 0, 7: 0}
+    
+
+
+# Check feasibility
+
+
+```python
+feasibility = mvc_feasibility(g, random_solution)
+print(f'The feasibility is {feasibility}\n')
+```
+
+    The feasibility is False
+    
+
+
+# Fix broken constraints
+
+
+```python
+fixed_solution = mvc_easy_fix(g, random_solution)
+print(f'After fixation, the solution is {fixed_solution}\n')
+```
+
+    After fixation, the solution is {0: 1, 1: 1, 2: 1, 3: 0, 4: 1, 5: 1, 6: 1, 7: 0}
+    
+
+
+
+```python
+feasibility = mvc_feasibility(g, fixed_solution)
+print(f'The feasibility is {feasibility}\n')
+```
+
+    The feasibility is True
+    
+
+
+# Visualisation
+
+
+```python
+fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(12.8,4.8))
+colour_b4fix = ['C1' if random_solution[n] else 'C0' for n in g.nodes]
+colour_fixed = ['C1' if fixed_solution[n] else 'C0' for n in g.nodes]
+nx.draw(g, pos, node_color=colour_b4fix, ax=ax0)
+ax0.set_title(f'Vertex cover in orange before fix')
+nx.draw(g, pos, node_color=colour_fixed, ax=ax1)
+ax1.set_title(f'Vertex cover in orange after fix')
+```
+
+
+
+
+    Text(0.5, 1.0, 'Vertex cover in orange after fix')
+
+
+
+
+    
+![png](README_files/README_22_1.png)
+    
+
+
+# Calculate energy
+
+
+```python
+energy = mvc_energy(g, fixed_solution)
+print(f'The energy is {energy}')
+```
+
+    The energy is 3.2102004750256556
+
diff --git a/examples/mvc/README_files/README_22_1.png b/examples/mvc/README_files/README_22_1.png
new file mode 100644
index 0000000000..deeb904926
Binary files /dev/null and b/examples/mvc/README_files/README_22_1.png differ
diff --git a/examples/mvc/README_files/README_7_0.png b/examples/mvc/README_files/README_7_0.png
new file mode 100644
index 0000000000..fd70caee3b
Binary files /dev/null and b/examples/mvc/README_files/README_7_0.png differ
diff --git a/examples/mvc/main.py b/examples/mvc/main.py
new file mode 100755
index 0000000000..7edb75a0a6
--- /dev/null
+++ b/examples/mvc/main.py
@@ -0,0 +1,73 @@
+"""Minimum Vertex Cover"""
+import csv
+import networkx as nx
+import argparse
+from qibo import callbacks, hamiltonians, models
+
+from mvc import qubo_mvc, qubo_mvc_penalty, mvc_feasibility, mvc_easy_fix, mvc_energy, hamiltonian_mvc
+from utils import binary2spin, spin2QiboHamiltonian
+
+parser = argparse.ArgumentParser()
+parser.add_argument('--filename', default='./mvc.csv', type=str)
+
+def load_csv(filename:str):
+    """ Load graph from csv file
+    """
+    with open(filename, 'r', newline='') as f:
+        reader = csv.reader(f)
+        data = [[int(row[0]), int(row[1]), float(row[2])] for row in reader]
+    
+    nodes = [k0 for k0, k1, v in data if k0 == k1]
+    edges = [[k0, k1] for k0, k1, v in data if k0 != k1]
+    
+    g = nx.Graph()
+    g.add_nodes_from(nodes)
+    g.add_edges_from(edges)
+    
+    node_weights = {k0: {'weight': v} for k0, k1, v in data if k0 == k1}
+    edge_weights = {(k0, k1): {'weight': v} for k0, k1, v in data if k0 != k1}
+
+    nx.set_edge_attributes(g, edge_weights)
+    nx.set_node_attributes(g, node_weights)
+
+    return g
+
+
+
+def main(filename: str):
+    
+    print(f'Load a graph from {filename} and make it a QUBO')
+    g = load_csv(filename)
+    penalty = qubo_mvc_penalty(g)
+    linear, quadratic = qubo_mvc(g, penalty=penalty)
+    
+    print('A random solution with seed 1234 must be infeasible')
+    import numpy as np
+    rng = np.random.default_rng(seed=1234)
+    random_solution = {i: rng.integers(2) for i in g.nodes}
+    feasibility = mvc_feasibility(g, random_solution)
+    assert not feasibility, 'The random solution should be infeasible.'
+    
+    print('Make a naive fix to the violation of the constraint')
+    fixed_solution = mvc_easy_fix(g, random_solution)
+    feasibility = mvc_feasibility(g, fixed_solution)
+    assert feasibility, 'The fixed solution should be feasible.'
+    
+    print('Calculate the energy of the solution')
+    energy = mvc_energy(g, fixed_solution)
+
+    print('Construct a hamiltonian based on the QUBO representation')
+    h, J, _ = binary2spin(linear, quadratic)
+    h = {k: -v for k, v in h.items()}
+    J = {k: -v for k, v in J.items()}
+    ham = spin2QiboHamiltonian(h, J)
+
+    print('Construct a hamiltonian directly from a networkx graph')
+    ham = hamiltonian_mvc(g)
+
+    print('done.')
+    
+if __name__ == "__main__":
+    # by defualt, test on the mvc.csv in the same directory
+    args = parser.parse_args()
+    main(args.filename)
diff --git a/examples/mvc/mvc.csv b/examples/mvc/mvc.csv
new file mode 100644
index 0000000000..b3f026c7cc
--- /dev/null
+++ b/examples/mvc/mvc.csv
@@ -0,0 +1,17 @@
+0,0,0.4696822179283675
+1,1,0.6917392308135916
+2,2,0.7130420958314817
+3,3,0.49342677332980056
+4,4,0.49661600571433673
+5,5,0.55601135361347
+6,6,0.2831095711244086
+7,7,0.18737823232636885
+0,1,0.98602725407379
+1,2,0.935853268681996
+2,3,0.23080434023289664
+4,5,0.25521731195623476
+5,6,0.37096743935296606
+6,7,0.04408120693490547
+0,2,0.5619060764177991
+4,7,0.48402095290884917
+1,6,0.7106584134580318
\ No newline at end of file
diff --git a/examples/mvc/mvc.py b/examples/mvc/mvc.py
new file mode 100644
index 0000000000..ca80b0006f
--- /dev/null
+++ b/examples/mvc/mvc.py
@@ -0,0 +1,138 @@
+import networkx as nx
+from utils import binary2spin, spin2QiboHamiltonian
+
+
+def qubo_mvc(g:nx.Graph, penalty=None):
+    """Given a graph g, return the QUBO of Minimum Vertex Cover (MVC) problem.
+    
+    If penalty is not given, it is inferred from g.
+
+    Args:
+        g (nx.Graph): graph
+        penalty (float): penalty weight of the constraints
+    
+    Returns:
+        linear (Dict[Any, float]): linear terms
+        quadratic (Dict[(Any, Any), float]): quadratic terms
+
+    """
+
+    q = {(k, k): v for k, v in g.nodes(data='weight')}
+    
+    if penalty is None:
+        penalty = qubo_mvc_penalty(g)
+
+    for s, d in g.edges:
+
+        q[(s, s)] -= penalty
+        q[(d, d)] -= penalty
+
+        if s > d:
+            s, d = d, s
+            
+        if (s, d) not in q:
+            q[(s, d)]  = penalty
+        else:
+            q[(s, d)] += penalty
+
+    linear    = {      k0: v for (k0, k1), v in q.items() if k0 == k1}
+    quadratic = {(k0, k1): v for (k0, k1), v in q.items() if k0 != k1 and v != 0}
+
+    return linear, quadratic
+
+
+def qubo_mvc_penalty(g:nx.Graph):
+    """Find appropriate penalty weight to ensure feasibility
+    
+    Args:
+        g (nx.Graph): graph
+    
+    Returns:
+        penalth (float): penalty
+
+    """
+
+    # For NISQ devices, you cannot randomly choose a positive penalty weight, otherwise
+    # you get infeasible solutions very likely. If all weight equals to 1,
+    # set penalty to the largest degree will ensure high feasibility.
+
+    highest_degree = max(sum(g.nodes[k]['weight'] for k in g[i].keys()) for i in g.nodes())
+    return highest_degree
+
+
+def mvc_feasibility(g:nx.Graph, solution):
+    """Check if the solution violates the constraints of the problem 
+    
+    Args:
+        g (nx.Graph): graph
+        solution (Dict[Any, int]): the solution
+    
+    Returns:
+        feasibility (bool): whether the solution meet the constraints
+    """
+    for k0, k1 in g.edges:
+        if solution[k0] == solution[k1] == 0:
+            return False
+    return True
+
+
+def mvc_energy(g:nx.Graph, solution):
+    """Calculate the energy of the solution on a MVC problem
+
+    The result is based on the assumption that soltuion is feasible.
+    
+    Args:
+        g (nx.Graph): graph
+        solution (Dict[Any, int]): the solution
+    
+    Returns:
+        energy (float): the energy of the solution
+    """
+    return sum(solution[k] * v for k, v in g.nodes.data('weight'))
+
+
+def mvc_easy_fix(g:nx.Graph, solution):
+    """Naively fix violation in an out-of--place manner
+    
+    Args:
+        g (nx.Graph): graph
+        solution (Dict[Any, int]): the solution
+    
+    Returns:
+        solution (Dict[Any, int]): fixed solution
+    """
+    
+    t = {k: v for k, v in solution.items()}
+
+    for k0, k1 in g.edges:
+        if t[k0] == t[k1] == 0:
+            t[k0] = 1
+            
+    return t
+
+
+def hamiltonian_mvc(g: nx.Graph, penalty: float=None, dense: bool=False):
+    """Given a graph g, return the hamiltonian of Minimum Vertex Cover (MVC)
+    problem.
+    
+    If penalty is not given, it is inferred from g.
+
+    Args:
+        g (nx.Graph): graph
+        penalty (float): penalty weight of the constraints
+        dense (bool): sparse or dense hamiltonian
+    
+    Returns:
+        ham: qibo hamiltonian
+
+    """
+
+    linear, quadratic = qubo_mvc(g, penalty)
+    
+    h, J, _ = binary2spin(linear, quadratic)
+    h = {k: -v for k, v in h.items()}
+    J = {k: -v for k, v in J.items()}
+    
+    ham = spin2QiboHamiltonian(h, J, dense=dense)
+
+    return ham
diff --git a/examples/mvc/utils.py b/examples/mvc/utils.py
new file mode 100644
index 0000000000..2acc638281
--- /dev/null
+++ b/examples/mvc/utils.py
@@ -0,0 +1,113 @@
+from typing import Dict, List, Tuple, Union, Any, Optional, Generator, Iterator, Hashable
+
+def binary2spin(linear: Dict[int, float],
+                quadratic: Dict[Tuple[int, int], float],
+                offset: float=0,):
+    """Convert binary model to spin model
+
+        Please remember to put a negative sign to h and J after using this
+        function , if you are going to form a mamiltonian from the spin
+        model. Hamiltonians usually have a leading negative sign in them,
+        but QUBOs don't.
+
+    Args:
+        linear (dict): linear term of the binary model
+        quadratic (dict): quadratic term of the binary model
+        offset (float): offset of the binary model
+
+    Returns:
+        h (dict): bias of the spin model
+        J (dict): interaction of the spin model
+        offset (float): offset of the spin model
+    """
+
+    h = dict((x, 0.5 * w) for x, w in linear.items())
+
+    J = []
+    for (x, y), w in quadratic.items():
+        J.append(((x, y), 0.25 * w))
+        h[x] += 0.25 * w
+        h[y] += 0.25 * w
+    J = dict(J)
+
+    offset += 0.5 * sum(linear.values())
+    offset += 0.25 * sum(quadratic.values())
+
+    return h, J, offset
+
+
+def spin2binary(h: Dict[int, float], 
+                J: Dict[Tuple[int, int], float],
+                offset: float=0,):
+    """Convert spin model to binary model
+
+        Please remember to put a negative sign to h and J before using this 
+        function if you extract them from a hamiltonian.  Hamiltonians
+        usually have a leading negative sign in them, but QUBOs don't.
+
+    Args:
+        h (dict): bias of the spin model
+        J (dict): interaction of the spin model
+        offset (float): offset of the spin model
+
+    Returns:
+        linear (dict): linear term of the binary model
+        quadratic (dict): quadratic term of the binary model
+        offset (float): offset of the binary model
+    """
+
+    linear = dict((s, 2. * bias) for s, bias in h.items())
+
+    quadratic = []
+    for (s, t), bias in J.items():
+        quadratic.append(((s, t), 4. * bias))
+        linear[s] -= 2. * bias
+        linear[t] -= 2. * bias
+    quadratic = dict(quadratic)
+
+    offset -= sum(linear.values())
+    offset += sum(quadratic.values())
+
+    return linear, quadratic, offset
+
+def spin2QiboHamiltonian(h: Dict[int, float], 
+                         J: Dict[Tuple[int, int], float],
+                         dense: bool=True,):
+    """Convert spin model to qibo Hamiltonian
+
+        Mixer is not included.
+
+        Please remember to put a negative sign to h and J if you get h and J
+        from binary2spin.  Hamiltonians usually have a leading negative sign
+        in them, but QUBOs don't.
+
+    .. math::
+        H = - \\sum_{i=0}^N \\sum _{j=0}^N J_{ij} Z_i Z_j - \\sum_{i=0}^N h_i Z_i.
+
+    Args:
+        h (dict): bias of the spin model
+        J (dict): interaction of the spin model
+        dense (bool): If ``True`` it creates the Hamiltonian as a
+            :class:`qibo.core.hamiltonians.Hamiltonian`, otherwise it creates
+            a :class:`qibo.core.hamiltonians.SymbolicHamiltonian`.
+
+    Returns:
+        linear (dict): linear term of the binary model
+        quadratic (dict): quadratic term of the binary model
+        offset (float): offset of the binary model
+    """
+
+    from qibo import hamiltonians
+    from qibo.symbols import Z
+
+    symbolic_ham  = sum(Z(k, commutative=True) * v for k, v in h.items())
+    symbolic_ham += sum(Z(k0, commutative=True) * Z(k1, commutative=True) * v
+                                                   for (k0, k1), v in J.items())
+    symbolic_ham = -symbolic_ham
+
+    ham = hamiltonians.SymbolicHamiltonian(symbolic_ham)
+
+    if dense:
+        return ham.dense
+    else:
+        return ham
diff --git a/examples/test_examples.py b/examples/test_examples.py
index 57aa0bc1ff..6de7cd366e 100644
--- a/examples/test_examples.py
+++ b/examples/test_examples.py
@@ -302,3 +302,11 @@ def test_grover_example3(nqubits, num_1):
     sys.path[-1] = path
     os.chdir(path)
     run_script(args, script_name="example3.py")
+
+
+def test_mvc():
+    args = locals()
+    path = os.path.join(base_dir, "mvc")
+    sys.path[-1] = path
+    os.chdir(path)
+    run_script(args)