AbstractRecursive equations for the efficient computation of nearly orthogonal, C0 continuous p‐version approximation functions (“max‐orthogonal”) are developed in one dimension and used to construct hexagonal (“brick”) element approximation functions. A modification to the functions is presented to simplify adaptive p‐refinement, resulting in “semi‐orthogonal” approximation functions. These functions are found to produce Laplace equation stiffness matrices, as well as mass matrices, with substantially lower condition numbers compared to standard Legendre‐based modal functions. Sample problems are presented using a conjugate gradient matrix solver to compare the matrix solution computation times using max‐orthogonal, semi‐orthogonal, and standard modal functions. Automatic adaptive polynomial refinement is also discussed and applied to these sample problems. The semi‐orthogonal functions are found to be an efficient and practical means of achieving very well‐conditioned p‐version finite element matrices.