syms U sigma epsilon E nu u(x,y) v(x,y) w(x,y) x y h

v=-h*diff(w,y);
u=-h*diff(w,x);

epsilon(1,1)=diff(u,x);
epsilon(2,2)=diff(v,y);
epsilon(1,2)=1/2*(diff(u,y)+diff(v,x));
epsilon(2,1)=epsilon(1,2);

sigma(1,1)=E/(1-nu^2)*(epsilon(1,1)+nu*epsilon(2,2));
sigma(2,2)=E/(1-nu^2)*(epsilon(2,2)+nu*epsilon(1,1));
sigma(1,2)=E/(1+nu)*epsilon(1,2);
sigma(2,1)=sigma(1,2);

U=E/(1-nu^2)*(1/2*epsilon(1,1)^2+1/2*epsilon(2,2)^2+nu*epsilon(1,1)*epsilon(2,2)+(1-nu)*epsilon(1,2)^2);
simplify(U)