syms a x1 x2 x3 real
u = a*[x1^2*x2; x3+2*x2^2; -x2*x3];
x = [x1; x2; x3];

G = zeros(3,3)*a;
for i=1:3
    for j=1:3
        G(i,j) = diff(u(i),x(j));
    end
end

a = 0.1;
% r = linspace(0,1,10);
% 
% X = zeros(3,length(r)^3);
% U = zeros(size(X));
% k = 0;
% for i1=1:length(r)
%     x1 = r(i1);
%     for i2=1:length(r)
%         x2 = r(i2);
%         for i3=1:length(r)
%             x3 = r(i3);
%             k = k+1;
%             X(:,k) = [x1 x2 x3]';
%             U(1,k) = eval(u(1));
%             U(2,k) = eval(u(2));
%             U(3,k) = eval(u(3));
%         end
%     end
% end

U = zeros(10,10);
V = zeros(10,10);
W = zeros(10,10);

for k=1:6
    switch k
        case 1
            X = linspace(0,1,10)'*ones(1,10);
            Y = X';
            Z = zeros(10,10);
        case 2
            X = linspace(0,1,10)'*ones(1,10);
            Y = X';
            Z = ones(10,10);
        case 3
            X = linspace(0,1,10)'*ones(1,10);
            Y = zeros(10,10);
            Z = X';
        case 4
            X = linspace(0,1,10)'*ones(1,10);
            Y = ones(10,10);
            Z = X';
        case 5
            X = zeros(10,10);
            Y = linspace(0,1,10)'*ones(1,10);
            Z = Y';
        case 6
            X = ones(10,10);
            Y = linspace(0,1,10)'*ones(1,10);
            Z = Y';

    end
    for i=1:10
        for j=1:10
            x1 = X(i,j);
            x2 = Y(i,j);
            x3 = Z(i,j);
            U(i,j) = eval(u(1));
            V(i,j) = eval(u(2));
            W(i,j) = eval(u(3));
        end
    end
    mesh(X+U,Y+V,Z+W)
    hold on
end
    
E = 1/2*(G + G');
W = 1/2*(G - G');

w = [W(3,2) W(1,3) W(2,1)]';
