Numerical computation of multiscale systems can be quite challenging, often requiring a very small discretization to adequately capture the fast component. In this talk I present a closure approximation method for the slow component of a two-timescale system of ODEs, replacing the coupled fast component with a simple parametrization. We use the fluctuation-dissipation relation, a result from statistical mechanics to predict perturbation response using observations of the unperturbed dynamics, to incorporate first-order response of the fast component in this parametrization. Results will be shown for the Lorenz 96 system for various coupling and forcing regimes.