The quadratic vibronic coupling model is an important computational tool for simulating photoelectron spectra involving strongly coupled electronic states in polyatomic molecules. However, recent work has indicated the need for higher order terms, with most of the initial studies focusing on molecules with symmetry-required degeneracies. In this study we report an extension of our approach for constructing fully quadratic representations of bound electronic states coupled by conical intersections, which allows for the inclusion of higher order terms, demonstrated here employing a quartic expansion. Procedures are developed that eliminate unphysical behavior for large displacements, a problem likely to be an endemic to anharmonic expansions. Following work on representing dissociative electronic states, Lagrange multipliers are used to constrain the constructed representation to reproduce exactly the energy, energy gradients, and/or derivative couplings at specific points, or nodes, in nuclear coordinate space. The approach is illustrated and systematically studied using the four lowest electronic states of triazolyl, (CH) 2 N 3 .