Accurate, reliable and fast numerical modeling methods are required to design the optimum radial refractive index profile for single and multimode fibers to give specific dispersion characteristics prior to or even obviating costly experimental work. Such profiles include graded index and multiple concentric cladding layers. In this paper, a new numerical method is introduced which enables the derivatives of the propagation coefficient to be calculated analytically up to the third order of a single mode or multimode weakly guiding optical fiber with an arbitrary radial refractive index profile. These quantities are required to determine the group delay, τg, chromatic dispersion, D, and dispersion slope of the fiber. The expansion of the modal fields in terms of Laguerre-Gauss polynomials in the Galerkin method offers certain benefits. In particular, due to simplicity of the basis functions it is possible to carry out further analytical work on the results such as repeated differentiation of the matrix equation resulting from the Galerkin method to define up to the third-order derivatives of the propagation coefficients with respect to wavelength. This avoids approximation errors inherent in numerical differentiation, giving better accuracy and, at the same time, significantly reduces the computation time. A computer program was developed to demonstrate the proposed method for single and multimode fibers with radially arbitrary refractive index profiles. The paper provides simulation results to validate the approach.