%eval(['[L,dL]=',func,'(x);'])
d=diag(.0001*abs(x));
d=diag(.0001*abs(x));
d(4,4)=.000001;
d(5,5)=.00001;
d(6,6)=.00001;
d(7,7)=.0000001;
d(8,8)=.0000001;
d(9,9)=.000001;
d(11,11)=1e-11;
n=length(x);
f=func;
%  H=zeros(n);
  for i=4;
    for j=1:i;
      eval(['H(i,j)=(',f,'(x+d(:,i)+d(:,j),par)-',f,'(x+d(:,i)-d(:,j),', ..
            'par)-',f,'(x-d(:,i)+d(:,j),par)+',f,'(x-d(:,i)-d(:,j),par))', ..
            '/(4*d(i,i)*d(j,j));']);
    end;
  end;
  H=H+tril(H,-1)';
