function d2 = secder(A,k,l)
% this function finds the second derivatives of lambda 
% for one element (at row k, column l) of the matrix a
% with respect to all other element. The function 
% returns a matrix of second derivatives: 
% d^2(lambda)/(daij,dakl)
n = size(A,1);
% find and sort the eigenvalues and vectors of A
[w,lambda]= eig(A);
v=inv(conj(w));
lambda = diag(lambda);
[lambda,index] = sort(lambda);
lambda = flipud(lambda);
index = flipud(index);
w=w(:,index);
v=v(index,:);
v=v.';
% calculate sensitivities
for i=1:n;
   senmat=conj(v(:,i))*(w(:,i).');
   svec(:,i) = senmat(:);
end
s1=svec(:,1);
for m = 2:n
   scalesens(:,m-1)=svec(:,m)/(lambda(1)-lambda(m));
end
% the line below is modified from Caswell 2001
if n>2 vecsum=sum(scalesens.'); else vecsum = scalesens.'; end;
%compute second derivatives
for i=1:n
   for j=1:n
      x1=(l-1)*n+1 + (i-1);
      x2=(j-1)*n+1 + (k-1);
      d2(i,j)= s1(x1)*vecsum(x2)+s1(x2)*vecsum(x1);
   end
end
d2=real(d2);