This paper concerns the numerical stability of an eikonal transformation based splitting method which is highly effective and efficient for the numerical solution of paraxial Helmholtz equation with a large wave number. Rigorous matrix analysis is conducted in investigations and the oscillation-free computational procedure is proven to be stable in an asymptotic sense. Simulated examples are given to illustrate the conclusion.