function yh=punishhessian(x,r)
syms x1 x2 x3;
y=1000-x1^2-2*x2^2-x3^2-x1*x2-x1*x3+r*((x1^2+x2^2+x3^2-25)^2+(8*x1+14*x2+7*x3-56)^2);
yh=[diff(y,x1,2) diff(diff(y,x1),x2) diff(diff(y,x1),x3);
    diff(diff(y,x2),x1) diff(y,x2,2) diff(diff(y,x2),x3);
    diff(diff(y,x3),x1) diff(diff(y,x3),x2) diff(y,x3,2)];
x1=x(1);
x2=x(2);
x3=x(3);
yh=subs(yh);