In this paper, a numerical method is presented to solve fractional boundary value problems. In fractional calculus, the modelling of natural phenomenons is best described by fractional differential equations. So, it is important to formulate efficient and accurate numerical techniques to solve fractional differential equations. In this article, first, we introduce $\psi$-shifted Chebyshev polynomials then project these polynomials to formulate $\psi$-shifted Chebyshev operational matrices. Finally, these operational matrices are used for the solution of fractional boundary value problems. The convergence is analysed. It is observed that solution of non-integer order differential equation converges to corresponding solution of integer order differential equation. Finally, the efficiency and applicability of method is tested by comparison of the method with some other existing methods.