A novel numerical technique based on the orthogonal polynomial approximation with minimum error-variance is presented for analyzing optical fibers with arbitrary refractive index profile. The developed method is verified to provide reliable solution in the case where the existing method based on the Chebyshev polynomial expansion shows instability. The usefulness of the proposed method lies in a number of desirable features viz., procedural simplicity, complete linearization of the problem, numerical stability, and more, that its application brings.