Describing how resource consumption rates depend on resource density, conventionally termed "functional responses," is crucial to understanding the population dynamics of trophically interacting organisms. Yet, accurately determining the functional response for any given pair of predator and prey remains a challenge. Moreover, functional responses are potentially complex surfaces in multidimensional space, where resource density is only one of several factors determining consumption rates. We explored how three sources of error can be addressed in the design and statistical analysis of functional response experiments: ill-chosen spacing of prey densities, heteroscedastic variance in consumption rates, and non-independence of parameters of the function describing prey consumption in relation to prey density and additional environmental factors. We generated extensive, virtual data sets that simulated feeding experiments in which both prey density and environmental temperature were varied, and for which the true, deterministic functional response surface was known and realistic variance had been added. We compared eight different methods of functional response fitting, one of which stood out as best performing. We subsequently tested several conclusions from the simulation study against experimental data of zooplankton feeding on algae across a broad range of temperatures. We summarize our main findings in three best practice guidelines for the experimental estimation of functional response surfaces, of which the second is the most important: (1) space prey densities logarithmically, starting from very low densities; (2) log-transform prey consumption data prior to fitting; and (3) fit a multivariate functional response surface to all data (including all prey densities and other factors, in our case temperature) in a single step. We also observed that functional response surfaces were fitted more accurately and precisely than their component parameters. The latter occurred because parameter estimates were non-independent, which is an inevitable feature of fitting complex nonlinear functions to data: A given response surface can often be described with near-equal accuracy by multiple parameter combinations. We therefore conclude that fitted functional response models perform better at optimizing the fit of the overall response surface than at determining how component parameters, such as the attack rate or handling time, depend on environmental factors such as temperature.