Consider the Hilbert space (H,〈• , •〉) equipped with the indefinite inner product[u,v]=v*J u,u,v∈ H, where J is an indefinite self-adjoint involution acting on H. The Krein space numerical range WJ(T) of an operator T acting on H is the set of all the values attained by the quadratic form [Tu,u], with u ∈H satisfying [u,u]=± 1. We develop, implement and test an alternative algorithm to compute WJ(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.