We propose a new family of functionally-fitted block θ-methods for numerically solving ordinary differential equations which integrates a chosen set of linearly independent functions exactly. The advantage of such variable coefficient methods is that the basis functions can be chosen to exploit specific properties of the problem that may be known in advance. The basic theory for the proposed methods is established. First, we derive a sufficient condition to ensure the existence of the functionally-fitted block θ-methods, and discuss the independence on integration time for a set of separable basis functions. We then obtain some basic characteristics of the methods by Taylor series expansions, and show that the order of accuracy of r-dimensional functionallyfitted block θ-method is at least r for ordinary differential equations. Numerical experiments are conducted to illustrate the efficiency of the functionally-fitted block θ-methods. 1. Introduction. Consider the following ordinary differential equations (ODEs) y (t) = f (t, y(t)), t ∈ [t 0 , T ], y(t 0) = y 0. (1) Suppose that f : [t 0 , T ] × R → R is continuous and satisfies a Lipschitz condition |f (t, y) − f (t, z)| ≤ L|y − z|, ∀t ∈ [t 0 , T ], y, z ∈ R such that there exists a unique solution. We shall assume that y(t) has continuous derivatives on [t 0 , T ] of any order needed. Numerous classical numerical methods including Runge-Kutta methods, multistep methods and block methods have