clc
clear all;
syms x y


S1 = [0, 200]';
S2 = [400, 400]';

h(1,1)=sqrt((x-S1(1))^2 + (y-S1(2))^2);
h(2,1)=sqrt((x-S2(1))^2 + (y-S2(2))^2);
h(3,1)=atan2((x-S1(2)),(y-S1(1)));
h(4,1)=atan2((x-S2(2)),(y-S2(1)));


f11 = symfun(h(1,1), [x y]);
f21 = symfun(h(2,1), [x y]);
f31 = symfun(h(3,1), [x y]);
f41 = symfun(h(4,1), [x y]);

value1=1;
value2=2;

D1=hessian(f11,[x y]);
D1_=D1(value1,value2);
D1_=[D1_(1,1) D1_(1,2) 0 0;...
    D1_(2,1) D1_(2,2) 0 0;...
    0 0 0 0;...
    0 0 0 0];

D2=hessian(f21,[x y]);
D2_=D2(value1,value2);
D2_=[D2_(1,1) D2_(1,2) 0 0;...
    D2_(2,1) D2_(2,2) 0 0;...
    0 0 0 0;...
    0 0 0 0];

D3=hessian(f31,[x y]);
D3_=D3(value1,value2);
D3_=[D3_(1,1) D3_(1,2) 0 0;...
    D3_(2,1) D3_(2,2) 0 0;...
    0 0 0 0;...
    0 0 0 0];

D4=hessian(f41,[x y]);
D4_=D4(value1,value2);
D4_=[D4_(1,1) D4_(1,2) 0 0;...
    D4_(2,1) D4_(2,2) 0 0;...
    0 0 0 0;...
    0 0 0 0];