function [ff,jac,hes,cub,df,djac,dhes]=hsf3(x,p)
d=p(1); h1=p(2); h2=p(3); theta=p(4); 
z=x(1); k=x(2); n1l=x(3); n1=x(4); n2=x(5); i=x(6); 
%
N=n1-(d/2)*(n1-n1l)^2;
f=z*k^theta*(h1*N^(1-theta) + h2*n2^(1-theta));
ff=f-i;
if nargout > 1;
  tem=[d*(n1-n1l);1-d*(n1-n1l)];
  jac=[f/z;
       theta*f/k;
       (1-theta)*z*k^theta*h1*N^(-theta)*tem;
       (1-theta)*z*k^theta*h2*n2^(-theta);
       -1];
end;
if nargout > 2;
  hes=zeros(6);
  hes(1,2:5)=jac(2:5)'/z;
  hes(2,2)=(theta-1)*jac(2)/k;
  hes(2,3:5)=theta*jac(3:5)'/k;
  hes(3:4,3:4)=(1-theta)*z*k^theta*h1*(-theta*N^(-theta-1)* ..
              tem*tem'-N^(-theta)*[d,-d;-d,d]);
  hes(5,5)=-theta*jac(5)/n2;
  hes=triu(hes)+triu(hes,1)';
end;

if nargout > 3;

  % cb = vec( dhes(x,p)/dx )

  cb=zeros(6);                       % wrt z
  cb(2:5,2:5)=hes(2:5,2:5)/z;
  cub(:,1)=cb(:);

  cb=zeros(6);                       % wrt k
  cb(1,2:5)=hes(2,2:5)/z;
  cb(2,2)=(theta-2)*hes(2,2)/k;
  cb(2,3:5)=(theta-1)*hes(2,3:5)/k;
  cb(3:4,3:4)=theta*hes(3:4,3:4)/k;
  cb(5,5)=theta*hes(5,5)/k;
  cb=triu(cb)+triu(cb,1)';
  cub(:,2)=cb(:);

  cb=zeros(6);                       % wrt n1l
  cb(1,2:5)=hes(3,2:5)/z;
  cb(2,2)=cub(9,2);
  cb(2,3:5)=cub(15:17,2)';
  cb(3,3)=(1-theta)*z*k^theta*h1*(-theta*(-theta-1)*N^(-theta-2)*tem(1)^3 ..
    +theta*N^(-theta-1)*2*d*d*(n1-n1l)-(-theta)*N^(-theta-1)*d*tem(1));
  cb(3,4)=(1-theta)*z*k^theta*h1*(-theta*(-theta-1)*N^(-theta-2)*tem(1)^2*  ..
    tem(2)+theta*N^(-theta-1)*d*(1-2*d*(n1-n1l))+ ..
    (-theta)*N^(-theta-1)*d*tem(1));
  cb(4,4)=(1-theta)*z*k^theta*h1*(-theta*(-theta-1)*N^(-theta-2)*tem(2)^2*  ..
    tem(1)-theta*N^(-theta-1)*2*(1-d*(n1-n1l))*d- ..
    (-theta)*N^(-theta-1)*d*tem(1));
  cb=triu(cb)+triu(cb,1)';
  cub(:,3)=cb(:);
  
  cb=zeros(6);                       % wrt n1
  cb(1,2:5)=hes(4,2:5)/z;
  cb(2,2)=cub(10,2);
  cb(2,3:5)=cub(21:23,2)';
  cb(3,3)=cub(16,3);
  cb(3,4)=cub(22,3);
  cb(4,4)=(1-theta)*z*k^theta*h1*(-theta*(-theta-1)*N^(-theta-2)*tem(2)^3 ..
    +theta*N^(-theta-1)*2*(1-d*(n1-n1l))*d -(-theta)*N^(-theta-1)*d*tem(2));
  cb=triu(cb)+triu(cb,1)';
  cub(:,4)=cb(:);
  
  cb=zeros(6);                       % wrt n2
  cb(1,2:5)=hes(5,2:5)/z;
  cb(2,2)=cub(11,2);
  cb(2,3:5)=cub(27:29,2)';
  cb(5,5)=(-theta-1)*hes(5,5)/n2;
  cb=triu(cb)+triu(cb,1)';
  cub(:,5)=cb(:);
  cub(:,6)=zeros(36,1);

end;
if nargout > 4;

  % df = df(x,p)/dp

  df=[z*k^theta*h1*(1-theta)*N^(-theta)*(-(n1-n1l)^2/2);
      z*k^theta*N^(1-theta);
      z*k^theta*n2^(1-theta);
      z*h1*(log(k)-log(N))*k^theta*N^(1-theta)+  ..
        z*h2*(log(k)-log(n2))*k^theta*n2^(1-theta)];  

end;

if nargout > 5;

  % djac = djac(x,p)/dp

  djac =[df'/z;
         theta*df'/k+[0,0,0,f/k];
         [jac(3)/d,jac(3)/h1,0,-jac(3)/(1-theta)+jac(3)*log(k/N)];
         [-jac(4)*(n1-n1l)/tem(2),jac(4)/h1,0, ..
             -jac(4)/(1-theta)+jac(4)*log(k/N)];
         [0,0,jac(5)/h2,-jac(5)/(1-theta)+jac(5)*log(k/n2)];
         zeros(1,4)]; 
end;

if nargout > 6;

  % dhes = vec( dhes(x,p)/dp )

  dh=zeros(6);
  dh(1,2:5)=djac(2:5,1)'/z;
  dh(2,2)=(theta-1)*djac(2,1)/k;
  dh(2,3:5)=theta*djac(3:5,1)'/k;
  dh(3:4,3:4)=hes(3,3)/d*[1,-1;-1,1];
  dh(5,5)=-theta*djac(5,1)/n2;
  dh=triu(dh)+triu(dh,1)';

  dhes(:,1)=dh(:);

  dh=zeros(6);
  dh(1,2:5)=djac(2:5,2)'/z;
  dh(2,2)=(theta-1)*djac(2,2)/k;
  dh(2,3:5)=theta*djac(3:5,2)'/k;
  dh(3:4,3:4)=hes(3:4,3:4)/h1;
  dh=triu(dh)+triu(dh,1)';

  dhes(:,2)=dh(:);

  dh=zeros(6);
  dh(1,2:5)=djac(2:5,3)'/z;
  dh(2,2)=(theta-1)*djac(2,3)/k;
  dh(2,3:5)=theta*djac(3:5,3)'/k;
  dh(5,5)=hes(5,5)/h2;
  dh=triu(dh)+triu(dh,1)';

  dhes(:,3)=dh(:);

  dh=zeros(6);
  dh(1,2:5)=djac(2:5,4)'/z;
  dh(2,2)=jac(2)/k+(theta-1)*djac(2,4)/k;
  dh(2,3:5)=jac(3:5)'/k+theta*djac(3:5,4)'/k;
  dh(3:4,3:4)=-hes(3:4,3:4)/(1-theta)-(1-theta)*z*(k/N)^theta*h1* ..
     ((theta/N*tem*tem'+[d,-d;-d,d])*log(k/N)+tem*tem'/N); 
  dh(5,5)=hes(5,5)/theta-theta*djac(5,4)/n2;
  dh=triu(dh)+triu(dh,1)';

  dhes(:,4)=dh(:);

end;
