We consider the numerical solution of the Bagley–Torvik equation Ay″(t) + BD3/2*y(t) + Cy(t) = f(t), as a prototype fractional differential equation with two derivatives. Approximate solutions have recently been proposed in the book and papers of Podlubny in which the solution obtained with approximate methods is compared to the exact solution. In this paper we consider the reformulation of the Bagley–Torvik equation as a system of fractional differential equations of order 1/2. This allows us to propose numerical methods for its solution which are consistent and stable and have arbitrarily high order. In this context we specifically look at fractional linear multistep methods and a predictor–corrector method of Adams type.