So this may be due to confusion about how mesh() works. Mesh produces a 2D surface in a 3D volume, with the Z-axis being the value of the matrix at each X/Y point. If your solution's are truly 3D (a value defined at each X/Y/Z point), mesh cannot really plot any more than one 2D slice at a time (at least cleanly). A good alternative would be to represent error as color in a 3D plot, using slice, but you won't be able to see every point at once.
You could also present the data a few other ways: 1. animated or subplotted 2D frames, plotting 1 2D "slice" at a time, 2. A single Mesh() showing the average or maximum error along the Z axis at each X/Y, and 3. Using something like histogram(err(:)) to see the breakdown of all the errors.
x = linspace( xstart, xend, nthroot(numel(zz),3) ); y = x;
and you can vectorize (compute all at once) the error generation by:
err = abs(u-v)
The if branch is useless, as error will always be >=0.
