nu = 0.3
G = 100;
del = (1-2*nu)/(2*(1-nu));
x = p(1,:);
y = p(2,:);
vx= 2*G*x.*y.^2.*(-1+y).^2.*(2*x-1).*(x-1);
vy= 2*G*x.^2.*y.*(2*y-1).*(-1+y).*(x-1).^2;
u1 = -vy + del*vx;
u2 = vx + del*vy;
U = [u1';u2'];
Umax = norm(U,inf)
Uhmax = norm(Uh,inf)
maxerror = norm(U-Uh,inf)