This paper considers the optimization of time-varying multivariable control systems with quadratic cost function which may have time-varying weighting matrices. The problem is first formulated in the context of functional analysis which leads to the optimal control being denned by a functional equation. An approximate solution of this equation is then obtained by replacing the kernel of the equation by a double Fourier series expansion. It is shown that the resulting sub-optimal control may be made to approach the absolute optimum as closely as desired by extending the double Fourier series, and bounds on the norm of the error in the control and the cost function are derived for any finite expansion. Application of the method to a particular system shows that Chebyshev basis functions lead to a better approximation than the classic sinusoidal functions. Finally, the use of the method as an adjunct to iterative optimization based on the contraction mapping algorithm is discussed briefly.