syms x y r theta nu A
nu = 0.3; A = 1;
r = sqrt(x^2+y^2);
theta = acos(x/r);
u1 = A*sqrt(r/2/pi)*(1-2*nu+sin(theta/2)^2)*cos(theta/2);
u2 = A*sqrt(r/2/pi)*(2-2*nu-cos(theta/2)^2)*sin(theta/2);
exx = simplify(diff(u1,x));
eyy = simplify(diff(u2,y));
exy = (simplify(diff(u1,y)+diff(u2,x)))/2;
pstrain = eig([exx,exy;exy,eyy])
fcontour(pstrain(1),[-0.2,0.2,0.1,0.5],'fill','on')

%

syms E11 E12 E22
Lagrange = [E11,E12;E12,E22];
eq1 = [1,0]*Lagrange*[1,0]'==(5/4-1)/2;
eq2 = [0,1]*Lagrange*[0,1]'==(5/4-1)/2;
eq3 = [-1,1]*Lagrange*[-1,1]' == (1-2)/4
[E11,E12,E22] = solve([eq1,eq2,eq3],[E11,E12,E22])

%

syms x1 x2 r e11 e22 e12 u1 u2 real
syms f(x2) g(x1)
r = sqrt(x1^2+x2^2)
e11 = x2^2/r^3;
e22 = x1^2/r^3;
e12 = - x1*x2/r^3;
simplify(diff(e11,x2,2) + diff(e22,x1,2) - 2*diff(diff(e12,x1),x2))
u1 = int(e11,x1) + f(x2)
u2 = int(e22,x2) + g(x1)
eq = diff(u1,x2) + diff(u2,x1) - 2*e12 ==0