For higher order elements, I refine each element a few times so I have more points to work with. If I just need to visualize the solution for myself.
Let's use quadratic Lagrange elements as examples. You need the mesh data, points
p and triangles
t, also the numerical solution $u_h$. For visualization purpose, we merely need the nodal values from the numerical solution using quadratic Lagrange elements (quadratic Lagrange has edge DoFs and node DoFs). Same for the higher order elements, we only need the nodal DoF's value to visualize the solution.
For example: your solution column vector is
u which is gonna be the $z$ value of the
surf plot. For quadratic element, the most natural data structure is to assign the first
NNodes = size(p,1) rows the value of nodal DoFs, then
NEdges = size(e,1) rows the value of edge DoFs. Because you use triangular mesh, we can use MATLAB's trisurf to do this:
h = trisurf(t, p(:,1), p(:,2), u(1:size(p,1))');
colormap's range can be manually assigned, see MATLAB's documents here. There are some built-in schemes like jet, you can use an $m\times 3$ RGB matrix as well (you can retrieve the color matrix by typing
c = colormap; then check what
c is like).
Here is another very helpful MATLAB doc on how to manually setting the
colormap for a patch object (
trisurf are both using patch to draw things): Coloring Mesh and Surface Plots
Here is what the numerical solution using quadratic elements is like for $-\Delta u = 1$ with homogeneous Dirichlet boundary data on an L-shape domain. There are no discontinuous gap between the nodes like your plot yields. Unless you use $hp$-Discontinuous Galerkin, your plot is wrong for usually continuous finite element solutions for any elliptic equations.