Reaction-diffusion equations have been widely used to describe biological pattern formation. Nonuniform steady states of reaction-diffusion models correspond to stationary spatial patterns supported by these models. Frequently these steady states are not unique, which correspond to various spatial patterns observed in biology. Traditionally, time-marching methods or steady state solvers based on Newton’s method were used to compute such solutions. However, the solution that any of these methods leads to highly depends on the initial condition/guess. In this talk, I present a systematic method to compute multiple nonuniform steady states for reaction-diffusion models and determine the dependence on model parameters. The method is based on homotopy continuation techniques and multigrid methods. We apply the method to two classic reaction-diffusion models and compare our results with available theoretical analysis in the literature. The first is the Schnakenberg model that has been used to describe biological pattern formation due to diffusion-driven instability. The second is the Gray-Scott model which was proposed in 1980’s to describe autocatalytic glycolysis reactions. In each case, our method uncovers many, if not all, nonuniform steady states and their stabilities. We also compared our computational results with analytical results in the literature and the comparison suggests some errors in prior results obtained using asymptotic analysis.