Multidisciplinary engineering systems are usually modeled by coupling software components that were developed for each discipline independently. The use of disparate solvers complicates the optimization of multidisciplinary systems and has been a long-standing motivation for optimization architectures that support modularity. The individual discipline feasible (IDF) formulation is particularly attractive in this respect. IDF achieves modularity by introducing optimization variables and constraints that effectively decouple the disciplinary solvers during each optimization iteration. Unfortunately, the number of variables and constraints can be significant, and the IDF constraint Jacobian required by most conventional optimization algorithms is prohibitively expensive to compute. Furthermore, limited-memory quasi-Newton approximations, commonly used for large-scale problems, exhibit linear convergence rates that can struggle with the large number of design variables introduced by the IDF formulation. In this work, we show that these challenges can be overcome using a reduced-space inexact-Newton-Krylov algorithm. The proposed algorithm avoids the need for the explicit constraint Jacobian and Hessian by using a Krylov iterative method to solve the Newton steps. The Krylov method requires matrix-vector products, which can be evaluated in a matrix-free manner using second-order adjoints. The Krylov method also needs to be preconditioned, and a key contribution of this work is a novel and effective preconditioner that is based on approximating a monolithic solution of the (linearized) multidisciplinary system. We demonstrate the efficacy of the algorithm by comparing it with the popular multidisciplinary feasible formulation on two test problems.