In a bounded domain, the solution of linear homogeneous time fractional parabolic equation is known to exhibit polynomial type decay rate (the so-called Mittag-Leffler stability) over time, which is quite different from the exponential decay of classical parabolic equation. We firstly use the finite element method or finite difference method to discretize the parabolic equation in space to obtain fractional ordinary differential equation, and then use fractional linear multistep method (F-LMM) to discretize in time to obtain a fully discretized schemes. We prove that the strongly $A$-stable F-LMM method combined with appropriate spatial discretization can accurately maintain the long-term optimal algebraic decay rate of the original continuous equation. Numerical examples are included to confirm the correctness of our theoretical analysis.