This paper investigates an interface-coupled fractional subdiffusion model, featuring two subdiffusion equations in adjacent domains connected by an interface allowing bidirectional energy transfer. The fractional derivative, accounting for long-term medium effects, introduces challenges in theoretical analysis and computational efficiency. We propose a partitioned time-stepping algorithm using higher-order extrapolations on the interface term to decouple the system with improved temporal accuracy, combined with finite element spatial approximations. Rigorous theoretical analysis demonstrates unconditional stability and optimal $L^2$ norm error estimates, supported by several numerical experiments.