function dirfld(xmin,xmax,ymin,ymax,nx,ny,f);
%
%   The direction field is plotted at (nx)x(ny) equally spaced
%   points in the area defined by (xmin,xmax)x(ymin,ymax).
%
%   Input parameters:
%         xmin,xmax   minimum and maximum values of x to use
%         ymin,ymax   minimum and maximum values of y to use
%         nx          number of points in the x-direction
%         ny          number of points in the y-direction
%         f           name of the function defining the ODE
%
%
%   The user must supply the function f(x,y) defining the
%   first order ODE dy/dx = f(x,y).  
%   The syntax is  
%            function yprime=f(x,y) 
%   where x and y are arrays of points and yprime is the corresponding
%   slope; i.e. yprime contains dy/dx at (x,y).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Written by Deborah Sulsky, Department of Mathematics and Statistics,
%     University of New Mexico, Albuquerque, NM 87131
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Clear figure window
    clf
%
%   Determine the spacing between points   
        hx=(xmax-xmin)/nx;
        hy=(ymax-ymin)/ny;
%   Make a grid of x and y values, using this spacing
        x=xmin+hx/2:hx:xmax-hx/2;
        y=ymin+hy/2:hy:ymax-hy/2;
        [X,Y]=meshgrid(x,y);
%   Make a grid of ones the same size as the xy-grid
        OO=ones(size(X));
%   Evaluate the slope of the direction field at the
%   xy-mesh points
        F=feval(f,X,Y);
%   Scale the line segments so they all have the same length
        scale=sqrt(1+F.*F);
        F=F./scale;
        OO=OO./scale;
%   Plot the direction field using quiver
        quiver(X,Y,OO,F,0.5,'.');
%   Set the plot frame
        axis([xmin,xmax,ymin,ymax]);
