Since the early 1970s, the study of Black–Scholes (BS) partial differential equations (PDEs) under the Efficient Market Hypothesis (EMH) has been a subject of active research in financial engineering. It has now become obvious, even to casual observers, that the classical BS models derived under the EMH framework fail to account for a number of realistic price evolutions in real-time market data. An alternative approach to the EMH framework is the Fractal Market Hypothesis (FMH), which proposes better and clearer explanations of market behaviours during unfavourable market conditions. The FMH involves non-local derivatives and integral operators, as well as fractional stochastic processes, which provide better tools for explaining the dynamics of evolving market anomalies, something that classical BS models may fail to explain. In this work, using the FMH, we derive a time-fractional Black–Scholes partial differential equation (tfBS-PDE) and then transform it into a heat equation, which allows for ease of implementing a high-order numerical scheme for solving it. Furthermore, the stability and convergence properties of the numerical scheme are discussed, and overall techniques are applied to pricing European put option problems.