In this dissertation, fitted operator and fitted mesh numerical methods are developed for singularly perturbed parabolic differential and differential-difference equations with boundary turning points. These problems are typically modeled by partial differential equations and differential-difference equations in which a small pa- rameter ε multiplies the highest-order derivative and the convection coefficient vanishes at the boundary, giv- ing rise to boundary turning points. Due to the presence of the perturbation parameter, the solution exhibits a multi-scale character, often featuring a thin transition layer where it varies abruptly, while away from the layer the solution behaves smoothly and changes gradually. Standard numerical methods on uniform meshes are generally unable to provide reliable approximations in such cases. The proposed methods in each chapter are ε-uniform, achieve better order of convergence and provide reliable solutions. The fitted operator numerical method is developed based on both the exponentially fitted operator approach and the non-standard finite dif- ference method on uniform meshes. For the fitted mesh schemes, the well-known piecewise-uniform Shishkin mesh is employed to accurately resolve the layers. Additionally, a post-processing technique (Richardson ex- trapolation) is used to enhance the accuracy of some of the developed methods. To validate the theoretical findings and demonstrate the accuracy of the proposed approaches, extensive numerical results are presented in tables and graphs for various mesh sizes, perturbation parameters ε, and orders of turning point p. The stability and convergence of the schemes are rigorously established using the discrete minimum principle, consistency arguments, and a priori error bounds, leading to ε-uniform convergence results. All simulations in this dissertation are performed using MATLAB. The results demonstrate the effectiveness of the proposed schemes and show improved accuracy and convergence compared to some existing methods.