Consider the Hilbert space (H; h ; i) equipped with the indefinite inner product OEu; v D v J u, u; v 2 H, where J is an indefinite self-adjoint involution acting on H. The Krein space numerical range W J .T / of an operator T acting on H is the set of all the values attained by the quadratic form OET u; u, with u 2 H satisfying OEu; u D˙1. We develop, implement and test an alternative algorithm to compute W J .T / in the finite dimensional case, constructing 2 by 2 matrix compressions of T and their easily determined elliptical and hyperbolical numerical ranges. The numerical results reported here indicate that this method is very efficient, since it is faster and more accurate than either of the existing algorithms. Further, it may yield easy solutions for the inverse indefinite numerical range problem. Our algorithm uses an idea of Marcus and Pesce from 1987 for generating Hilbert space numerical ranges of matrices of size n.