In this paper we develop a numerical approach to a fractional-order differential Linear Complementarity Problem (LCP) arising in pricing European and American options under a geometric Lévy process. The LCP is first approximated by a nonlinear penalty fractional Black-Scholes (fBS) equation. We then propose a finite difference scheme for the penalty fBS equation. We show that both the continuous and the discretized fBS equations are uniquely solvable and establish the convergence of the numerical solution to the viscosity solution of the penalty fBS equation by proving the consistency, stability and monotonicity of the numerical scheme. We also show that the discretization has the 2nd-order truncation error in both the spatial and time mesh sizes. Numerical results are presented to demonstrate the accuracy and usefulness of the numerical method for pricing both European and American options under the geometric Lévy process.