  A simple MATLAB implementation of hierarchical shape functions on 2D
rectangles is explained and available for download. Global shape functions are
ordered for a given polynomial degree according to the indices of the nodes,
edges, or elements to which they belong. For a uniform p-refinement, the
hierarchical structure enables an effective assembly of mass and stiffness
matrices. A solution of a boundary value problem is approximated for various
levels of uniform h and p refinements.
