The current documentation provides the code for the tau polynomials as
# Tau polynomials
tau_basis = xbasis.derivative_basis(2)
p1 = dist.Field(bases=tau_basis)
p2 = dist.Field(bases=tau_basis)
p1['c'][-1] = 1
p2['c'][-2] = 2
There is no explanation of the choice of tau_basis or the manual setting of the elements of p1 and p2. This should be clarified. I can do the work and open a pull request with some pointers on the explanation.