When　the　fluctuation　of　option　price　is　regarded　as　a　fractal　transmission　system　and　the　stock　price　follows　a　Lévy　distribution,　a　time-space　fractional　option　pricing　model　(TSFOPM)　is　obtained.　Then　we　discuss　the　numerical　simulation　of　the　TSFOPM.　A　discrete　implicit　numerical　scheme　with　a　second-order　accuracy　in　space　and　a　2−γ　order　accuracy　in　time　is　constructed,　where　γ　is　a　transmission　exponent.　The　stability　and　convergence　of　the　obtained　numerical　scheme　are　analyzed.　Moreover,　a　fast　bi-conjugate　gradient　stabilized　method　is　proposed　to　solve　the　numerical　scheme　in　order　to　reduce　the　storage　space　and　computational　cost.　Then　a　numerical　example　with　exact　solution　is　presented　to　demonstrate　the　accuracy　and　effectiveness　of　the　proposed　numerical　method.　Finally,　the　TSFOPM　and　the　above　numerical　technique　are　applied　to　price　European　call　option.　The　characteristics　of　the　fractional　option　pricing　model　are　analyzed　in　comparison　with　the　classical　Black–Scholes　(B-S)　model.　©　2018　Elsevier　Inc.