function ffs=DeriveDeterminant(C,Cs)
% DERIVEDETERMINANT Newton correction for Determinant root finding
%   ffs=DeriveDeterminant(C,Cs) computes for the two given matrices
%   C(lambda) and its derivative C'(lambda), lambda being a real
%   parameter, the Newton correction ffs=f/fs for the determinant
%   function f(C(lambda))=0.

n=length(C); lfs=0; 
for i=1:n
  [cmax,kmax]=max(abs(C(i:n,i)));
  if cmax==0                                    %  Matrix singular
    error('Matrix is singular')
  end
  kmax=kmax+i-1;
  if kmax~=i
    h=C(i,:); C(i,:)=C(kmax,:); C(kmax,:)=h;
    h=Cs(kmax,:); Cs(kmax,:)=Cs(i,:); Cs(i,:)=h;
  end
  lfs=lfs+Cs(i,i)/C(i,i);
  Cs(i+1:n,i)=(Cs(i+1:n,i)*C(i,i)-Cs(i,i)*C(i+1:n,i))/C(i,i)^2;
  C(i+1:n,i)=C(i+1:n,i)/C(i,i);
  Cs(i+1:n,i+1:n)=Cs(i+1:n,i+1:n)-Cs(i+1:n,i)*C(i,i+1:n)- ...
                  C(i+1:n,i)*Cs(i,i+1:n);
  C(i+1:n,i+1:n)=C(i+1:n,i+1:n)-C(i+1:n,i)*C(i,i+1:n);
end
ffs=1/lfs;