diff --git a/abp/graphstate.py b/abp/graphstate.py
index 31fc891..baa1c9d 100755
--- a/abp/graphstate.py
+++ b/abp/graphstate.py
@@ -49,7 +49,7 @@ class GraphState(object):
         self._add_node(self, *args, **kwargs)
 
     def _del_node(self, node):
-        """ Remove a node. TODO: this is a hack right now! """
+        """ Remove a node. TODO: this is a hack right now. """
         if not node in self.node:
             return
         del self.node[node]
@@ -57,6 +57,10 @@ class GraphState(object):
             del self.adj[k][node]
         del self.adj[node]
 
+    def del_qubit(self, node):
+        """ Remove a qubit. TODO: this is a hack right now. """
+        self._del_node(node)
+
     def _add_node(self, node, **kwargs):
         """ Add a node. By default, nodes are initialized with ``vop=``:math:`I`, i.e. they are in the :math:`|+\\rangle` state.
 
@@ -424,8 +428,7 @@ class GraphState(object):
         """ Represent as a string for quick debugging """
         s = ""
         for key in sorted(self.node.keys()):
-            s += "{}:  {}\t".format(
-                key, clifford.get_name(self.node[key]["vop"]).replace("YC", "-"))
+            s += "{}:  {}\t".format(key, clifford.get_name(self.node[key]["vop"]))
             if self.adj[key]:
                 s += str(tuple(self.adj[key].keys())).replace(" ", "")
             else:
diff --git a/abp/qi.py b/abp/qi.py
index 30db417..4d760f5 100644
--- a/abp/qi.py
+++ b/abp/qi.py
@@ -123,10 +123,18 @@ class CircuitModel(object):
         b = normalize_global_phase(other.state)
         return np.allclose(a, b)
 
+    def __setitem__(self, key, value):
+        """ Set a matrix element """
+        self.state[key] = value
+
+    def __getitem__(self, key):
+        """ Get a matrix element """
+        return self.state[key]
+
     def __str__(self):
         s = ""
         for i in range(self.d):
-            label = bin(i)[2:].rjust(self.nqubits, "0")
+            label = bin(i)[2:].rjust(self.nqubits, "0")[::-1]
             if abs(self.state[i, 0]) > 0.00001:
                 term = self.state[i, 0]
                 real_sign = " " if term.real>=0 else "-"
diff --git a/abp/static/index.html b/abp/static/index.html
index b9ceaa3..36c66ed 100644
--- a/abp/static/index.html
+++ b/abp/static/index.html
@@ -25,7 +25,7 @@
 <div id=node_info class=hidden> nothing </div>
 <div id=server_info class=hidden> </div>
 
-<div id=version>Version 0.4.20 develop mode</div>
+<div id=version>Version 0.4.20</div>
 
 <div id=node_data class=hidden>
 <div id=node_name></div>
diff --git a/doc/index.rst b/doc/index.rst
index b6fc530..6652b2e 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -45,32 +45,84 @@ If you want to modify and test ``abp`` without having to re-install, switch into
 Quickstart
 ----------------------------
 
-It's pretty easy to build a graph state, act some gates, and do measurements::
+Let's make a new ``GraphState`` object with a register of three qubits:
 
     >>> from abp import GraphState
-    >>> g = GraphState(5)
-    >>> for i in range(5):
-    ...     g.act_hadamard(i)
-    ... 
-    >>> for i in range(4):
-    ...     g.act_cz(i, i+1)
-    ... 
-    >>> print g 
-    0:  IA	(1,)
-    1:  IA	(0,2)
-    2:  IA	(1,3)
-    3:  IA	(2,4)
-    4:  IA	(3,)
-    >>> g.measure(2, "px")
-    0
+    >>> g = GraphState(3)
+
+All the qubits are initialized by default in the :math:`|+\rangle` state:
+
+    >>> print g.to_state_vector()
+    |000❭: 	 √1/8	+ i √0
+    |100❭: 	 √1/8	+ i √0
+    |010❭: 	 √1/8	+ i √0
+    |110❭: 	 √1/8	+ i √0
+    |001❭: 	 √1/8	+ i √0
+    |101❭: 	 √1/8	+ i √0
+    |011❭: 	 √1/8	+ i √0
+    |111❭: 	 √1/8	+ i √0
+
+We can also check the stabilizer tableau:
+
+    >>> print g.to_stabilizer()
+     0  1  2
+    ---------
+     X      
+        X   
+           X
+
+Or look directly at the vertex operators and neighbour lists:
+
+    >>> print g
+    0:  IA	-
+    1:  IA	-
+    2:  IA	-
+
+This representation might be unfamiliar. Each row shows the index of the qubit, then the _vertex operator_, then a list of neighbouring qubits. To understand vertex operators, read the original paper by Anders and Briegel.
+
+Let's act a Hadamard gate on the zeroth qubit -- this will evolve qubit ``0`` to the :math:`H|+\rangle = |1\rangle` state:
+
+    >>> g.act_hadamard(0)
+    >>> print g.to_state_vector()
+    |000❭: 	 √1/4	+ i √0
+    |010❭: 	 √1/4	+ i √0
+    |001❭: 	 √1/4	+ i √0
+    |011❭: 	 √1/4	+ i √0
+    
     >>> print g
-    0:  IA	(3,)
-    1:  ZC	(3,)
+    0:  YC	-
+    1:  IA	-
     2:  IA	-
-    3:  ZA	(0,1,4)
-    4:  IA	(3,)
 
-Working with GraphStates
+And now run some CZ gates:
+
+    >>> g.act_cz(0,1)
+    >>> g.act_cz(1,2)
+    >>> print g
+    0:  YC	-
+    1:  IA	(2,)
+    2:  IA	(1,)
+
+    >>> print g.to_state_vector()
+    |000❭: 	 √1/4	+ i √0
+    |010❭: 	 √1/4	+ i √0
+    |001❭: 	 √1/4	+ i √0
+    |011❭: 	-√1/4	+ i √0
+
+Tidy up a bit:
+
+    >>> g.del_node(0)
+    >>> g.act_hadamard(0)
+    >>> print g.to_state_vector()
+    |00❭: 	 √1/2	+ i √0
+    |11❭: 	 √1/2	+ i √0
+
+Cool, we made a Bell state. Incidentally, those those state vectors and stabilizers are real objects with methods, not just string-like representations of the state:
+
+
+
+
+GraphState API
 -------------------------
 
 The ``abp.GraphState`` class is the main interface to ``abp``.  
diff --git a/tests/test_qi.py b/tests/test_qi.py
index 421c09e..9add17b 100644
--- a/tests/test_qi.py
+++ b/tests/test_qi.py
@@ -169,3 +169,10 @@ def test_sqrt_notation(n=2):
     g = GraphState(range(n))
     g.act_circuit(c)
 
+def test_indexint():
+    """ Test that we can index into state vectors """
+    psi = qi.CircuitModel(0)
+    assert psi[0] == 1+0j
+    psi[0] = 42
+    assert psi[0] == 42
+