Stationary perturbation theory is concerned with finding the changes in the discrete energy levels and the changes in the corresponding energy eigenfunctions of a system, when the Hamiltonian of a system is changed by a small amount.  Let H = H₀ + H' = H₀ + λW.  H₀ is the unperturbed Hamiltonian whose eigenvalues E₀^p and eigenstates |Φ^p_i> are known.  Let {|Φ^p_i>} denote an orthonormal eigenbasis of H₀, H₀|Φ^p_i>= E₀^p|Φ^p_i>.  Here i denotes the degeneracy.  Assume that the matrix elements of H' in the eigenbasis of H₀ are small compared to the matrix elements of H₀.
<Φ^p_i|H'|Φ^p_i>  <<  <Φ^p_i|H₀|Φ^p_i> = E₀^p.
We write H' = λW with λ << 1 and <Φ^p_i|H₀|Φ^p_i> on the order of E₀^p.
We are looking for the eigenvalues E(λ) and the eigenstates |ψ(λ)> of H(λ) = H₀ + λW.
H(λ)|ψ_p> = E^p(λ)|ψ_p> or, keeping the dependence on λ in mind without specifically writing it down, H|ψ_p> = E^p|ψ_p>.
Since λW is small, we assume that E and |ψ> can be expanded as a power series in λ.

E^p = E₀^p + λE₁^p + λ²E₂^p + ... ,
|ψ_p> = |ψ_p⁰> + λ|ψ_p¹> + λ²|ψ_p²> + ... .
We may then write
(H₀ + λW)(|ψ_p⁰> + λ|ψ_p¹> + λ²|ψ_p²> + ...)
= (E₀^p + λE₁^p + λ²E₂^p + ...)(|ψ_p⁰> + λ|ψ_p¹> + λ²|ψ_p²> + ...).

This equation is must be valid over a continuous range of λ.  Therefore we equate coefficients of equal powers of λ on both sides to obtain a series of equations that represent successively higher orders of the perturbation.

• (H₀ - E₀^p)|ψ_p⁰> = 0
• (H₀ - E₀^p)|ψ_p¹> = (E₁^p - W)|ψ_p⁰>
• (H₀ - E₀^p)|ψ_p²> = (E₁^p - W)|ψ_p¹> + E₂^p|ψ_p⁰>
• (H₀ - E₀^p)|ψ_p³> = (E₁^p - W)|ψ_p²> + E₂^p|ψ_p¹> + E₃^p|ψ_p⁰>
• ...

(H₀ - E₀^p)|ψ_p⁰> = 0 implies that |ψ_p⁰> is a linear combination of unperturbed eigenfuctionns |Φ^p_i> with the corresponding eigenvalue E₀.  We choose <ψ_p⁰|ψ_p⁰> = 1.  |ψ_p^s>is not uniquely defined.  We can add an arbitrary multiple of |ψ_p⁰> to each |ψ_p^s> without affecting the left hand side of the above equations.  Most often this multiple is chosen so that <ψ_p⁰|ψ_p^s> = 0.  The perturbed ket is then not normalized.
We then have
<ψ_p⁰|(H₀ - E₀^p)|ψ_p^s> = <ψ_p⁰|(E₁^p - W)|ψ_p^s-1> + <ψ_p⁰|E₂^p|ψ_p^s-2> + ... + <ψ_p⁰|E_s^p|ψ_p⁰>,
or
0 = -<ψ_p⁰|W|ψ_p^s-1> + E_s^p,
or
E_s^p = <ψ_p⁰|W|ψ_p^s-1>.

To calculate the energy to sth order, we only need to know the state vector to order s - 1.

Perturbation theory for non-degenerate levels

First-order energy corrections

Consider a particular non-degenerate eigenvalue E₀ⁿ of H₀.  H₀|Φ_n> = E₀ⁿ|Φ_n>.  The other eigenvalues of H may or may not be degenerate.  We have |ψ_n⁰> = |Φ_n> and E₁ⁿ = <Φ_n|W|Φ_n>.
The first-order energy correction therefore is λE₁ⁿ = <Φ_n|H'|Φ_n>.
We have Eⁿ = E₀ⁿ + <Φ_n|H'|Φ_n> + O(λ²).

First-order eigenvector corrections

(H₀ - E₀^p)|ψ_p¹> = (E₁^p - W)|ψ_p⁰> = (E₁^p - W)|Φ_p>.  |ψ_p⁰> is an eigenstate of the unperturbed Hamiltonian.  We may expand
|ψ_p¹> = ∑_p'≠p,i b_p'ⁱ|Φ_p'ⁱ> in terms of the basis vectors |Φ_p'ⁱ>. 
In the expansion b_p = 0 because <ψ_p⁰|ψ_pⁱ> = 0.
∑_p'≠p,i(H₀ - E₀^p) b_p'ⁱ|Φ_p'ⁱ> = ∑_p'≠p,i(E₀^p' - E₀^p) b_p'ⁱ|Φ_p'ⁱ> = (E₁^p - W)|Φ_p>.

Multiply from the left by <Φ_p''ⁱ| we obtain
(E₀^p' - E₀^p) b_p''ⁱ = -<Φ_p''ⁱ|W|Φ_p>,
or
b_p'ⁱ = <Φ_p'ⁱ|W|Φ_p>/(E₀^p - E₀^p'),  p' ≠ p.
Therefore

|ψ_p> = |Φ_p> + ∑_p'≠p,i <Φ_p'ⁱ|H'|Φ_p>/(E₀^p - E₀^p')]|Φ_p'ⁱ> + O(λ²).

Second-order energy corrections

Since we have found the expression for the state vector to first order, we can now find the expression for the energy to second order.
E₂^p = <ψ_p⁰|W|ψ_p¹> = ∑_p'≠p,i b_p'ⁱ<Φ_p|W|Φ_p'ⁱ> = ∑_p'≠p,i <Φ_p'ⁱ|W|Φ_p><Φ_p|W|Φ_p'ⁱ>/(E₀^p - E₀^p').

E₂^p = ∑_p'≠p,i |<Φ_p'ⁱ|W|Φ_p>|²/(E₀^p - E₀^p').

E^p = E₀^p + <Φ_n|H'|Φ_n> + ∑_p'≠p,i |<Φ_p'ⁱ|W|Φ_p>|²/(E₀^p - E₀^p') + O(λ³)
= E₀^p + λE₁^p + λ²E₂^p + ... .
H = H₀ + H' = H₀ + λW.  In practice, after having derived the perturbation expansion, we often set λ = 1 and let H' = W be small.

Perturbation theory for degenerate levels

First-order energy corrections

Let H = H₀ + H' = H₀ + λW.  Let {|Φ^p_i>} denote an orthonormal eigenbasis of H₀, H₀|Φ_pⁱ> = E₀^p|Φ_pⁱ>.  Here i denotes the degeneracy.  Consider a particular degenerate eigenvalue E₀^p of H₀.  Assume that this eigenvalue is g-fold degenerate, i = 1,2,...,g.  To find E₁^p we use

(H₀ - E₀^p)|ψ_p¹> =  (E₁^p - W)|ψ_p⁰>, |ψ_p⁰> = ∑_ia_pⁱ|Φ_pⁱ >.

Multiplying from the left by <Φ_p^j| we obtain
0 = ∑_ia_pⁱ<Φ_p^j|(E₁^p - W)|Φ_pⁱ> = E₁^pa_p^j - ∑_ia_pⁱ<Φ_p^j|W|Φ_pⁱ >,
or
∑_i<Φ_p^j|W - E₁^pδ_ij|Φ_pⁱ>a_pⁱ = 0.

This is an eigenvalue equation for the operator W in the subspace E(0,p) of vectors with the eigenvalue of H₀ equal to E₀^p.  We find g eigenvalues E₁^p,i which may or may not be degenerate.  We then find the corresponding g eigenvectors.  To solve for the eigenvalues we set det(W - E₁^p) = 0 in the subspace E(0,p).  If we find a non-degenerate eigenvalue E₁^p,i, then the corresponding eigenvector is uniquely defined.  If the eigenvalue E₁^p,i is degenerate then the corresponding eigenvector is still not uniquely defined.  The degeneracy may or may not be removed in higher order.