Neutral particles play an important role on the performance of tokamak plasmas. In this dissertation, the original TEP methodology has been extended to take into account linearly (DP_1) and quadratically (DP_2) anisotropic distributions of angular fluxes for calculations of transmission probabilities. Three approaches, subdivision of optically thick regions, expansion of collision sources and the diffusion approximation, have been developed and implemented to correct effects of the preferential probability of collided neutrals escaping back across the incident surface. Solving the diffusion equation via the finite element method has been shown to be the most computationally efficient and accurate for a broader range of D/l by comparisons with Monte Carlo simulations. The average neutral energy (ANE) approximation has been developed and implemented into the GTNEUT code. The average neutral energy approximation has been demonstrated to be more accurate than the original local ion temperature (LIT) approximation for optically thin regions. The simulations of the upgraded GTNEUT code excellently agree with the DEGAS predictions in DIII-D L-mode and H-mode discharges, and the results of both the codes are in a good agreement with the experimental measurements.