Shooting schemes
for the shallow water equations
on a β-plane

Functions for solving the ODE

soln[ω_,k_,LR_,Lf_,L_]:=NDSolve[{
			v''[y]+(-(1+y/Lf)^2-k^2+ω^2-(k/ω)*(LR/Lf))v[y]==0,
		v[0]==0,
			v'[0]==1},
		v,{y,0,L}]

vAtL[ω_,k_,LR_,Lf_,L_]:=v[L]/.soln[ω,k,LR,Lf,L][[1]]

Some parameters we don't change much

{k,LR,Lf,L}={.1,1,10,1};

Try it!

trysoln=soln[-.001,k,LR,Lf,L]

ParametricPlot[{v[y]/. trysoln[[1]],y},{y,0,1},AspectRatio->2,Ticks->{{},Automatic}]

Try shooting for v=0 at y=L:

We see there are some high frequency solutions:

Plot[vAtL[ω,k,LR,Lf,L],{ω,-8,-0.001`}]

...but also there are some low frequency, quasi-geostrophic  solutions:

Plot[vAtL[ω,k,LR,Lf,L],{ω,-0.001`,-0.0001`}]

Let's find a root forω exactly:

aroot=FindRoot[vAtL[ω,k,LR,Lf,L],{ω,-.001,-.0008}]

The above does not work in Mathematica 8.0, so make a Interpolation of a Table:

desperate=Table[{ω,vAtL[ω,k,LR,Lf,L]},{ω,-.001,-.0008,.00001}]

aroot=FindRoot[Interpolation[desperate][ω],{ω,-.001,-.0008}]

Now find v[y] that satisfies the boundary conditions:

hitsoln=soln[ω/.aroot[[1]],k,LR,Lf,L]

ParametricPlot[{v[y]/. hitsoln[[1]],y},{y,0,1},AspectRatio->2,Ticks->{{},Automatic}]

Clearly the fun is just beginning, but the rest is up to you!