function [UU,VV]=dirfield(ff,x,y);
% DIRFIELD  plots a direction field
%
%  Example:  Assume you want the direction field of
%    dy
%    --  = x + 2*sin(y)      for -1<=x<=1, -3<=y<=3
%    dx
%
%  Write a function file "ff.m" as follows:
%     function dy = ff(x,y);
%     dy= x + 2*sin(y);
%
%  Then run:
%  >> dirfield('ff',-1:.2:1,-3:.5:3);
%
%  (The list "-1:.2:1" is the same as [-1 -.8 -.6 ....8 1].)

[XX,YY]=meshgrid(x,y);
ANG=atan(feval(ff,XX,YY));

UU=cos(ANG);
VV=sin(ANG);

quiver(XX,YY,UU,VV,.25,'.');
hold on;
quiver(XX,YY,-UU,-VV,.25,'.');
hold off;