In [14], the authors developed a new approach to the computation of the Hausdorff dimension of the invariant set of an iterated function system or IFS. In this paper, we extend this approach to incorporate high order approximation methods. We again rely on the fact that we can associate to the IFS a parametrized family of positive, linear, Perron-Frobenius operators $L_s$, an idea known in varying degrees of generality for many years. Although $L_s$ is not compact in the setting we consider, it possesses a strictly positive $C^m$ eigenfunction $v_s$ with eigenvalue $R(L_s)$ for arbitrary $m$ and all other points $z$ in the spectrum of $L_s$ satisfy $|z| le b$ for some constant $b < R(L_s)$. Under appropriate assumptions on the IFS, the Hausdorff dimension of the invariant set of the IFS is the value $s=s_*$ for which $R(L_s) =1$. This eigenvalue problem is then approximated by a collocation method at the extended Chebyshev points of each subinterval using continuous piecewise polynomials of arbitrary degree $r$. Using an extension of the Perron theory of positive matrices to matrices that map a cone $K$ to its interior and explicit a priori bounds on the derivatives of the strictly positive eigenfunction $v_s$, we give rigorous upper and lower bounds for the Hausdorff dimension $s_*$, and these bounds converge rapidly to $s_*$ as the mesh size decreases and/or the polynomial degree increases.