function xdot = f( s, x )

% Munk sound speed profile

eps = 0.00737;
c0  = 1500;
z   = x( 2 );
xt  = 2 * ( z - 1300 ) / 1300;  
c   = c0 * ( 1 + eps * ( xt - 1 + exp( -xt ) ) );
c2  = c^2;

% we also need derivatives of sound speed

dxtdz = 2 / 1300;
cz= c0 * eps * dxtdz * ( 1 - exp( -xt ) ); 
cr = 0;

% here's the RHS

xdot = zeros( 4, 1 );

xdot( 1 ) = c * x( 3 );
xdot( 2 ) = c * x( 4 );
xdot( 3 ) = -cr / c2;
xdot( 4 ) = -cz / c2;