In this paper, we describe in detail a scheme for the construction of highly accurate numerical solutions to Einstein's field equations in five and six spacetime dimensions corresponding to non-uniform black strings. The scheme consists of a sophistically adapted multi-domain pseudo-spectral method which incorporates a detailed understanding of the solution's behavior at the domain boundaries and at critical points... In particular, the five-dimensional case is exceedingly demanding as logarithmic terms appear which need to be treated with special care. The scheme resolves these issues and permits the investigation of unprecedentedly strong deformations of the black string horizon. As a consequence, we are able to study in detail the critical regime in phase diagrams displaying characteristic thermodynamic quantities such as mass and entropy. Our results show typical spiral curves in such diagrams which provides a strong support of previous numerical works. (read more)