We use the Constraint Satisfaction Problem (CSP) framework to model and solve the problem of designing substitution functions for substitution- permutation (SP) networks as proposed by Shannon for the architecture of ciphers. Many ciphers are designed using the SP pattern, and differ mainly by two parametrized functions: substitution and permutation. The most difficult of the two is the substitution function, which has to be nonlinear (a requirement that was difficult to define and quantify). Over time, researchers such as Nyberg, Pieprzyk and Matsui have proposed various metrics of nonlinearity that make the function robust to modern attacks. Before us, people have attempted various ways to design functions that respect these metrics. In the past people hand-picked substitution tables (S-boxes) by trying various values. Recently they use difficult to analyze constructs (such as Bent functions, spectral inversion, inverses in Galois fields) whose outputs are tested for nonlinearity. While efficient, such techniques are neither exhaustive (optimal), nor did they manage to generate better substitutions than the ones hand-picked in the past. We show that Matsui's nonlinearity requirement can be naturally modelled using CSPs. Based on a combination of existing CSP techniques and some new filtering operators that we designed specially for the new types of constraints, we manage to obtain better S-boxes than any previously published ones. The simplicity of the CSP framework and availability of general CSP solvers like ours, makes it easy for more people to design their own ciphers with easy to understand security parameters. Here we report on this new application of CSPs.