Option　pricing　models　are　often　used　to　describe　the　dynamic　characteristics　of　prices　in　financial　markets.　Unlike　the　classical　Black–Scholes　(BS)　model,　the　finite　moment　log　stable　(FMLS)　model　can　explain　large　movements　of　prices　during　small　time　steps.　In　the　FMLS,　the　second-order　spatial　derivative　of　the　BS　model　is　replaced　by　a　fractional　operator　of　order　α　which　generates　an　α-stable　Lévy　process.　In　this　paper,　we　consider　the　finite　difference　method　to　approximate　the　FMLS　model.　We　present　two　numerical　schemes　for　this　approximation:　the　implicit　numerical　scheme　and　the　Crank–Nicolson　scheme.　We　carry　out　convergence　and　stability　analyses　for　the　proposed　schemes.　Since　the　fractional　operator　routinely　generates　dense　matrices　which　often　require　high　computational　cost　and　storage　memory,　we　explore　three　methods　for　solving　the　approximation　schemes:　the　Gaussian　elimination　method,　the　bi-conjugate　gradient　stabilized　method　(Bi-CGSTAB)　and　the　fast　Bi-CGSTAB　(FBi-CGSTAB)　in　order　to　compare　the　cost　of　calculations.　Finally,　two　numerical　examples　with　exact　solutions　are　presented　where　we　also　use　extrapolation　techniques　to　achieve　higher-order　convergence.　The　results　suggest　that　the　proposed　schemes　are　unconditionally　stable　and　convergent,　and　the　FMLS　model　is　useful　for　pricing　options.　©　2020　Wiley　Periodicals,　Inc.