function unitary()
A = [1 3+i;3-i 4]
[X,D] = eig(A)
Y = X'
B = X*X'
X*D*X'-A
Y1 = [-5/(3-i);1]
z1 = A*Y1
z1 - D(1,1)*Y1
Y2 = [2/(3-i);1]
z2 = A*Y2
z2-D(2,2)*Y2
C=[(-i-3)/2; 1]
D=[(i-3)/2 1]
D*C
C'*C
Q=[(3+i)/5;1]
S=[(3-i)/5 1]
S*Q
Q'*Q
7/5