  In Evans function computations of the spectra of asymptotically
constant-coefficient linear operators, a basic issue is the efficient and
numerically stable computation of subspaces evolving according to the
associated eigenvalue ODE. For small systems, a fast, shooting algorithm may be
obtained by representing subspaces as single exterior products
\cite{AS,Br.1,Br.2,BrZ,BDG}. For large systems, however, the dimension of the
exterior-product space quickly becomes prohibitive, growing as $\binom{n}{k}$,
where $n$ is the dimension of the system written as a first-order ODE and $k$
(typically $\sim n/2$) is the dimension of the subspace. We resolve this
difficulty by the introduction of a simple polar coordinate algorithm
representing ``pure'' (monomial) products as scalar multiples of orthonormal
bases, for which the angular equation is a numerically optimized version of the
continuous orthogonalization method of Drury--Davey \cite{Da,Dr} and the radial
equation is evaluable by quadrature. Notably, the polar-coordinate method
preserves the important property of analyticity with respect to parameters.
