Let the smallest non-trivial solution of Pell equation, x^2-Dy^2=1, be denoted by (x_1, y_1). The Pell equations for D were systematically classified into several types with respect to the form of the polynomial relations (PR's) among (D, x_1, y_1). The key strategies for this analysis are the value of y_1 and form of the continued fraction expansion of √<D>. Among the solutions of Pell equation with D below 100 only four D's were found to have no other D below ten thousand connected through a PR. All the PR's were shown to be derived from a pair of the "master equations". The results obtained in this paper show an effectiveness of the proposed strategies for the systematic analysis of the chaotic behavior of the solutions of Pell equation