In this thesis, we study slowly rotating relativistic stars by modeling them as perfect fluids. The Hartle formalism is used to perturb a spherically symmetric static solution into slow rotation and the corresponding field equations to second order are presented. Five of the seven equations form a subsystem that decouples and can thus be solved independently. It is solely this system we focus on in this work. Two new variables are introduced to rewrite the field equations such that they become algebraic in one of the variables, the other one is used as a generating function and can be chosen to generate desired perfect fluid configurations. To adapt for numerical integration, the field equations are further reformulated as a coupled system of ordinary differential equations, which is solved numerically using a fourth order Runge-Kutta method. Two choices of the generating function are examined, namely a simple polynomial case as well as the Kuch2 III solution, and the related field equations are solved to find the perturbative functions. By studying physical quantities such as the pressure, energy density and speed of sound to zeroth order, we show that both of the studied configurations correspond to physically realistic solutions within a range of parameter values. Additionally, the shape of the zero pressure surface to second order is analyzed and found to be oblate for both cases, which coincides with the expected shape of a rotating fluid. We also investigate whether the Petrov type D condition is fulfilled, and the results indicate that both solutions fail to satisfy the condition and are thus not of Petrov type D. Moreover, it is shown that both the polynomial case and the perturbed Kuch2 III solution can be matched to a suitable asymptotically flat vacuum exterior, using the Darmois-Israel procedure. It should be noted that even though we are solely examining the slowly rotating limit, our results are in fact valid even up to fairly high rotational velocities. This indicates that the perturbative generating method of this thesis may be used to produce models describing physically reasonable configurations, e.g. rotating neutron stars.