function dxdt = v(x)
  % return velocity vector v(x) for cylinder flow at position x
  V0 = 1;
  a  = 1;

  r = sqrt(x(1)^2 + x(2)^2);
  theta = atan2(x(2),x(1));

  dxdt = [V0*(1 - a^2/r^2 * cos(2*theta));
         -V0*a^2/r^2 * sin(2*theta)];
	 
end