Load the package

Get["Eigen`"]

If you wish to use parallelization, load Eigen on each parallel kernal

ParallelNeeds["Eigen`"](*Unless your differential equations are very complicated, parallelization probably will not give a significant speed up*)

Harmonic Oscillator

The syntax of EigenNDSolve is very similar to that of NDSolve

?EigenNDSolve

Solve for the eigenvalues of the quantum harmonic oscillator on the domain (-10,10)

esys=EigenNDSolve[{-(1/2)ξ''[x]+(1/2)(x^2)ξ[x]==ω ξ[x],ξ[-10]==0,ξ[10]==0},ξ,{x,-10,10},ω];

EigenNDSolve returns the list {eigenvalues, eigenfunctions of (Subscript[u, 1]), eigenfunctionsof(Subscript[u, 2]),…}

(*In this case there is only one independent variable*)
Shallow[esys]

Conveniently plot the eigenspectrum using PlotSpectrum

?PlotSpectrum

PlotSpectrum[esys,{{0,20},{-0.8,0.6}}]
(*Click on the eigenvalues to display the eigenfunctions!*)

Any sort of homogeneous boundary conditions can be used

(*Neuman BCs*)
esysN=EigenNDSolve[{-(1/2)ξ''[x]+(1/2)(x^2)ξ[x]==ω ξ[x],ξ'[-10]==0,ξ'[10]==0},ξ,{x,-10,10},ω];
(*Periodic BCs*)
esysP=EigenNDSolve[{-(1/2)ξ''[x]+(1/2)(x^2)ξ[x]==ω ξ[x],ξ[-10]==ξ[10],ξ'[-10]==ξ'[10]},ξ,{x,-10,10},ω];

There is only limited checking of the supplied boundary conditions (since sometimes it is not obvious how many are needed) and output will be flawed if problem is over/under-specified.

PlotSpectrum[
EigenNDSolve[{-(1/2)ξ''[x]+(1/2)(x^2)ξ[x]==ω ξ[x],ξ[-10]==0},ξ,{x,-10,10},ω],
{{0,20},{-0.8,0.6}}]

One dimensional fluid eigenproblems

The Orr-Somerfeld equations arise from linearizing and Fourier analyzing a Navier-Stokes equilibrium with background flow U(x). α and β are the y and z wavenumbers respectively.

OrrSomer={(1/ReU)((α^4) u[x]+2 (α^2) (β^2) u[x]+(β^4) u[x]-2 (α^2) u''[x]-2 (β^2) u''[x]+Derivative[4][u][x])+I (α^3) u[x] U[x]+I α (β^2) u[x] U[x]+I ω (u''[x]-(α^2)  u[x]- (β^2) u[x])-I α U[x] u''[x]+I α u[x] U''[x]==0,(1/ReU)(ζ''[x]-(α^2) ζ[x]-(β^2) ζ[x])+I ω ζ[x]-I α U[x] ζ[x]-I β u[x] U'[x]==0};
BCs={u[-1]==0,u[1]==0,u'[-1]==0,u'[1]==0,ζ[-1]==0,ζ[1]==0};

Pouiseille flow U(x)=1-(x^2) is stable at low Reynolds number but becomes unstable above Re≈6000

PouiseilleEqns=Join[OrrSomer,BCs]/.{U->(x|->1-(x^2))}/.{α->1,β->0};
esysLR=EigenNDSolve[PouiseilleEqns/.ReU->1000,{u,ζ},{x,-1,1},ω];
esysHR=EigenNDSolve[PouiseilleEqns/.ReU->6000,{u,ζ},{x,-1,1},ω];
PlotSpectrum[esysLR,{{0,1},{-1,0.1}}]
PlotSpectrum[esysHR,{{0,1},{-1,0.1}}]

There are a couple of numerical effects to watch out for that are observable in the above figures. 
The first is the spurious eigenvalues, the first and second eigenvalues in both esysLR and esysLR are very large and not physical.
The second is the splitting of the lower eigenvalue branch in the high Reynolds number case. This can be improved by using more polynomials.

esysHR=EigenNDSolve[PouiseilleEqns/.ReU->6000,{u,ζ},{x,-1,1},ω,Npoly->130];
PlotSpectrum[esysHR,{{0,1},{-1,0.1}}]

Note that it is sometimes necessary to use higher than numerical precision to obtain satisfactory results (e.g., precision->30). This can be very slow!
More information on these methods and numerical challenges can be found in Dongarra et al. “Chebyshev tau QZ Algorithm Methods for Calculating Spectra of Hydrodynamic Stability Problems,” 1997.

The method used to construct operator matrices can be fairly slow for very complicated differential equations. If you find this to be a problem, setting useParallel->True can speed this up substantially. Note that the native Mathematica function Eigensystem is already parallelized for numerical precision calculations and in most cases this step is the most computationally intensive.