Cylindrical dielectric waveguides such as the optical fiber and photonic crystal fiber are very important passive devices in optical communication systems. There are many kinds of commercial software and methods of simulation at present. In this thesis, we proposed the following four methods to analyze arbitrarily profiled cylindrical dielectric waveguides: The first two methods are modified from published work while the last two methods are entirely developed by ourselves. 1. Cylindrical ABCD matrix method: We take the four continuous electromagnetic field components as main variables and derive the exact four-by-four matrix (with Bessel functions) to relate the four field vector within each homogeneous layer. The electromagnetic field components of the inner and outer layer can propagate toward one of the selected interface of our choice by using the method of ABCD matrix. We can then solve for the £]-value of the waveguide mode with this nonlinear inhomogeneous matrix equation. 2. Runge-Kutta method: Runge-Kutta method is mostly used to solve the initial value problems of the differential equations. In this thesis, we introduce the Runge-Kutta method to solve the first-order four-by-four nonlinear differential equation of the electromagnetic field components and find the £]-value of the cylindrical dielectric waveguides in a similar way depicted in method one. 3. Coupled Ez and Hz method: It uses the axial electromagnetic filed components to solve cylindrical dielectric waveguides. The formulation is similar to cylindrical ABCD matrix method, but it requires less variables then cylindrical ABCD matrix method. The numerical solution obtained from this method is most stable, but it is more complicated to derive harder to write the program. 4. Simple basis expansion method: The simple trigonometric functions (sine or cosine) are chosen as the bases of the horizontal coupled magnetic field equation derived from the second-order differential equation of the transverse magnetic field components. We do not select the horizontal coupling electric field because the normal component of the electric field is discontinuous on the interface. But the normal and tangential components of the magnetic field are continuous across the interfaces. The modal solution problem is converted to a linear matrix eigenvalue-eigenvector equation which is solved by the standard linear algebra routines. We will compare these four numerical methods with one another. The characteristics and advantage as well as the disadvantage of each method will be studied and compared in detail.