For a series of Markov processes we prove stochastic duality relations with duality functions given by orthogonal polynomials. This means that expectations with respect to the original process (which evolves the variable of the orthogonal polynomial) can be studied via expectations with respect to the dual process (which evolves the index of the polynomial). The set of processes include interacting particle systems, such as the exclusion process, the inclusion process and independent random walkers, as well as interacting diffusions and redistribution models of Kipnis-Marchioro-Presutti type. Duality functions are given in terms of classical orthogonal polynomials, both of discrete and continuous variable, and the measure in the orthogonality relation coincides with the process stationary measure.