This paper presents a new theoretical and computational framework for computing solutions of right classes for laminated composites using 2D p-version least squares finite element formulation incorporating the correct physics of interlamina behavior. At the interface between two laminas of dissimilar materials we have continuity of displacements u, v, stresses σyy, τxy, and strain εxx, while the stress σxx and the strains εyy and γxy are discontinuous. Thus, a finite element formulation, incorporating the physics of laminate behavior, would require interpolation of u, v, εxx, σyy and τxy instead of u, v, σxx, σyy and τxy which is generally the case in most mixed formulations. In the p-version LSFEF presented here, we interpolate u, v and σyy, τxy (εxx = ∂u/∂x is used to eliminate εxx as a variable) using appropriate p-version interpolations which would ensure correct interlamina behavior of these components. When the mating lamina properties are different, interlamina discontinuity of σxx, εyy and γxy is automatically generated due to dissimilar material properties of the laminas. In this formulation interlamina jumps in σxx, εyy and γxy do not constitute singularities, hence mesh refinements and higher p-levels are not needed in the vicinity of inter-lamina boundaries. The major thrust of this paper is to construct interpolations for the dependant variables that are of right classes in appropriate spaces so that a sequence of converged solutions in these spaces may be computed which, when converged, would yield a numerical solution that has exactly the same characteristics (in terms of continuity and differentiability) as the analytical or theoretical (strong) solution.