Global spectral methods offer the potential to compute solutions of partial differential equations numerically to very high accuracy. In this work, we develop a novel global spectral method for linear partial differential equations on cubes by extending the ideas of Chebop2 (Townsend, A. & Olver, S. (2015) The automatic solution of partial differential equations using a global spectral method. J. Comput. Phys., 299, 106-123) to the three-dimensional setting utilizing expansions in tensorized polynomial bases. Solving the discretized partial differential equation involves a linear system that can be recast as a linear tensor equation. Under suitable additional assumptions, the structure of these equations admits an efficient solution via the blocked recursive solver (Chen, M. & Kressner, D. (2020) Recursive blocked algorithms for linear systems with Kronecker product structure. Numer. Algorithms, 84, 1199-1216). In the general case, when these assumptions are not satisfied, this solver is used as a preconditioner to speed up computations.