The numerical simulation method for accurately solving the coupled generalized nonlinear Schrödinger equations (C-GNLSE) is essential for describing the dynamic behavior of ultrashort pulse propagating in optical fiber and developing the corresponding nonlinear fiber-optic devices. C-GNLSE in the normal picture is first mapped into the interaction picture by the representation transformation, and then, the two coupled nonlinear partial differential equations in the interaction picture are solved in frequency domain, with synchronous data updating in each iteration step, by using the vector form of Hult’s fourth-order Runge-Kutta iterative scheme. The proposed vector form algorithm of fourth-order Runge-Kutta in interaction picture (V-JH-RK4IP) is verified by using it in simulating the supercontinuum generation in high birefringence photonic crystal fiber. Moreover, the V-JH-RK4IP algorithm also exhibits the highest accuracy and computational efficiency as compared to other classical algorithms.