The aim of this thesis is to propose and analyze some efficient numerical methods for several kind of fractional differential equations. Some backgrounds for fractional differential equations will be first introduced in Chapter 1. In Chapter 2, compact finite difference schemes for the modified anomalous fractional sub-diffusion equation and fractional diffusion-wave equation are studied. Based on the idea of weighted and shifted Grunwald difference operator, we establish schemes ¨ with temporal and spatial accuracy order equal to two and four respectively. Numerical examples are carried out to support the theoretical results. In the next part, a high order finite difference scheme for one dimensional nonlinear fractional Klein-Gordon equations is first derived. The solvability of the difference system is discussed by the Leray-Schauder fixed point theorem, while the stability and L∞ convergence of the finite difference scheme are proved by the energy method. A two dimensional case subject to Neumann boundary conditions is also considered. The difficulty induced by the nonlinear term and the Neumann conditions is carefully handled in the proposed scheme. The stability and convergence of the finite difference scheme are also analyzed using the matrix form of the scheme. Numerical examples are provided to demonstrate the theoretical results of both cases. In Chapter 4, a high order compact exponential alternating direction implicit (ADI) scheme for two dimensional fractional convection-diffusion equations is proposed with O(τ 3−γ + h 4 1 + h 4 2 ) accuracy, where τ, h1, h2 are the temporal and spatial step sizes respectively. The convergence of the finite difference scheme is studied using its matrix form by the energy method. Difficulty arising from the convection term is overcome by analysis based on the eigenvalue decomposition of nonsymmetric tridiagonal matrices. Numerical examples are given to justify the theoretical results. In Chapter 5, a high order compact scheme for time fractional Fokker-Planck equations with variable convection is constructed. The scheme is studied using its matrix form by the energy method. We find that the difficulty arising from the variable coeffi- cient can be overcome by simple modifications of the coefficient matrices. The scheme is shown to be stable and convergent with order τ 2−α + h 4 which is higher than some recently studied schemes. Numerical examples are also provided. This thesis ends with some brief conclusions and future work.