You can solve this problem using mixed linear integer prrogramming, as follows:

- loop over the number n of colors
- for each such n, add n binary variables to each vertex and to each edge: bv[v,c] and be[e,c], where v is a vertex, e is an edge, and 0<=c<=n-1 is an integer. The meaning is: bv[v,c] is 1 iff the vertex v has color c, be[e,c] is 1 iff the edge e has color c.
- add the following constraints:
- for each vertex v, the sum of bv[v,c] for c<=n-1 is equal to 1 (meaning each vertex get exactely one color)
- for each edge, the sum of be[e,c] for c<=n-1 is equal to 1 (meaning each edge get exactly one color)
- for each vertex v, for each c<=n-1, bv[v,c] plus the sum of be[e,c] where e is an edge adjacent to v is at most 1 (coloring condition around vertices)
- for each edge e = (v,w), for each c<=n-1, bv[v,c]+bv[w,c]+be[e,c] is at most 1 (coloring condition around edges)

- try to solve, if you succeed, you get your coloring, if you get an exception, increase n by 1 and iterate

**EDIT** Here is a possible implementation;

```
sage: from sage.numerical.mip import MIPSolverException
sage: def total_chromatic_number(G, certificate=False):
....: nmax = len(G) + len(G.edges()) # trivial upper bound on the number of colors.
....: for n in range(1,nmax+1):
....: p = MixedIntegerLinearProgram()
....: bv = p.new_variable(binary=True)
....: be = p.new_variable(binary=True)
....: for v in G.vertices():
....: p.add_constraint(sum(bv[v,c] for c in range(n)) == 1)
....: for e in G.edges(labels=False):
....: p.add_constraint(sum(be[e,c] for c in range(n)) == 1)
....: for v in G.vertices():
....: for c in range(n):
....: p.add_constraint(bv[v,c] + sum(be[e,c] for e in G.edges_incident(v, labels=False)) <= 1)
....: for v,w in G.edges(labels=False):
....: for c in range(n):
....: p.add_constraint(bv[v,c] + bv[w,c] + be[(v,w),c] <= 1)
....: try:
....: p.solve()
....: if certificate:
....: bv_sol = p.get_values(bv)
....: be_sol = p.get_values(be)
....: coloration = {}
....: for v in G.vertices():
....: for c in range(n):
....: if bv_sol[v,c] == 1:
....: coloration[v] = c
....: for e in G.edges(labels=False):
....: for c in range(n):
....: if be_sol[e,c] == 1:
....: coloration[e] = c
....: return coloration
....: else:
....: return n
....: except MIPSolverException:
....: pass
```

We can test in on the Petersen graph, and see that its total chromatic number is 4:

```
sage: G = graphs.PetersenGraph()
sage: total_chromatic_number(G)
4
```

And get an admissible total coloration:

```
sage: total_chromatic_number(G, certificate=True)
{0: 0,
1: 1,
2: 3,
3: 0,
4: 1,
5: 3,
6: 0,
7: 0,
8: 2,
9: 2,
(0, 1): 3,
(0, 4): 2,
(0, 5): 1,
(1, 2): 0,
(1, 6): 2,
(2, 3): 2,
(2, 7): 1,
(3, 4): 3,
(3, 8): 1,
(4, 9): 0,
(5, 7): 2,
(5, 8): 0,
(6, 8): 3,
(6, 9): 1,
(7, 9): 3}
```