While general quantum many-body systems require exponential resources to be simulated on a classical computer, systems of non-interacting fermions can be simulated exactly using polynomially scaling resources. Such systems may be of interest in their own right, but also occur as effective models in numerical methods for interacting systems, such as Hartree-Fock, density functional theory, and many others. Often it is desirable to solve for equilibirum as well as non-equilibrium properties of inhomogeneous systems of many thousand constituent particles, rendering these simulations computationally costly despite their polynomial scaling. We demonstrate how this scaling can be improved by adapting methods based on matrix product states (MPS), which have been enormously successful for low-dimensional interacting quantum systems, to the case of free fermions. Compared to MPS simulations of interacting systems, our methods achieve an exponential speedup in the entanglement entropy of the state. We demonstrate their use to solve systems of up to one million sites with an effective MPS bond dimension of 10^15.