function [v,B3]=polyf(M)
% [v,B3]=polyf(M)
% Start with matrix M,
% Produce vector of coefficients v
% (negated and excluding leading)
% and 3-d array of (all) adjugate
% coefficients
     test=size(M);
     if (test(1)~=test(2)) 
         error('only square matrices are accepted here.');
     end%if
     n=test(1);
     B=eye(n);
     v=zeros(1,0);
     B3=zeros([test,0]);
% empty arrays used as basis
     for k=1:n
	B3=cat(3,B3,B);
        A=B*M;
        c=trace(A)/k;
	v=[v,c];
	B=A-c*eye(n);
      end%for
% last array should be zero,
% echoed for error detection
      B
% Results should be suppressed
end%function
% Rutgers teaching code by R. T. Bumby