We study computational theory and methods for finding multiple unstable solutions (corresponding to saddle points) to three types of nonlinear variational elliptic systems: cooperative, noncooperative, and Hamiltonian. We first propose a new Lorthogonal selection in a product Hilbert space so that a solution manifold can be defined. Then, we establish, respectively, a local characterization for saddle points of finite Morse index and of infinite Morse index. Based on these characterizations, two methods, called the local min-orthogonal method and the local min-max-orthogonal method, are developed and applied to solve those three types of elliptic systems for multiple solutions. Under suitable assumptions, a subsequence convergence result is established for each method. Numerical experiments for different types of model problems are carried out, showing that both methods are very reliable and efficient in computing coexisting saddle points or saddle points of infinite Morse index. We also analyze the instability of saddle points in both single and product Hilbert spaces. In particular, we establish several estimates of the Morse index of both coexisting and non-coexisting saddle points via the local min-orthogonal method developed and propose a local instability index to measure the local instability of both degenerate and nondegenerate saddle points. Finally, we suggest two extensions of an L-orthogonal selection for future research so that multiple solutions to more general elliptic systems such as nonvariational elliptic systems may also be found in a stable way.