This article is concerned with the problem of predicting a deterministic response function yo over a multidimensional domain T, given values of yo and all of its first derivatives at a set of design sites (points) in T. The intended application is to computer experiments in which yo is an output from a computer model of a physical system and each point in T represents a particular configuration of the input parameters. It is assumed that the first derivatives are already available (e.g., from a sensitivity analysis) or can be produced by the code that implements the model. A Bayesian approach in which the random function that represents prior uncertainty about yo is taken to be a stationary Gaussian stochastic process is used. The calculations needed to update the prior given observations of yo and its first derivatives at the design sites are given and are illustrated in a small example. The issue of experimental design is also discussed, in particular the criterion of maximizing the reduction in entropy, which leads to a kind of D optimality. It is shown that, for certain classes of correlation functions in which the intersite correlations are very weak, D-optimal designs necessarily maximize the minimim distance between design sites. A simulated annealing algorithm is described for constructing such maximin distance designs. An example is given based on a demonstration model of eight inputs and one output, in which predictions based on a maximin design, a Latin hypercube design, and two compromise designs are evaluated and compared.