A Lagrangian formulation for the valence bands of bulk zinc blende and wurtzite compound semiconductors provides a direct approach to determining derivative operator ordering at layer interfaces in the multiband k·P description of electronic states, in the envelope function approximation. The current continuity condition is obtained through a gauge-variation on the Lagrangian. The principle of least action together with the discretization of the action integral naturally leads into a finite element approach for the modeling and simulation of the multiband Schrödinger equations. Being a variational method, the wavefunctions and energy eigenvalues can be systematically improved with quadratic convergence. By including the Poisson Lagrangian, a self-consistent treatment of the Schrödinger-Poisson band-bending in arbitrarily doped structures is obtained. Numerical examples are presented for the valence band electronic states for quantum wells and superlattices. We have also coupled the finite element method with layer parameter optimizations to design heterostructures for lasers emitting at given wavelengths. It is in this sense that wavefunction engineering goes beyond the notion of bandgap engineering in providing the means for optimizing wavefunction localization in layers, selecting layer thicknesses, enhancing optical matrix elements, and laser gain.