function ay = fy(vx,vy)
  m = 1;
  b = 1;
  g = 10;
  [r,s] = size(vx);
  vter = sqrt(m*g/b);
  V = sqrt(vx.^2+vy.^2);
  ay = g*ones(r,s) - b/m*V.*vy;
endfunction