
function [] =  dfield(xylimits,npts,filename,scale)
%  File dfield.m plots derivative values as a field of slope vectors
%  for a 1st order ODE.
%  Follows approach by Acheson (1997), 'From Calculus to Chaos',
%  Oxford Press Ch. 3
%
%  xylimits = [xmin xmax ymin ymax]
%           = area over which derivative field is plotted
%  npts = number of grid points in x and y for plotting the slope field .
%  filename is a character string containing the name of the function file
%  used to compute the derivative (NEEDS .*, .^, etc., to operate on matrices).
%  scale adjusts length of slope vectors (typically 0.7).
%  (Vector lengths are scaled to be a factor=scale of the grid point spacing;
%  aspect ratio of plot not taken into account, but is okay.)
%
%
     xmin = xylimits(1);
     xmax = xylimits(2);
     ymin = xylimits(3);
     ymax = xylimits(4);
%
     dx = (xmax-xmin)/npts;     % spacing of grid points
     dy = (ymax-ymin)/npts;
%
[xm,ym] = meshgrid(xmin:dx:xmax,ymin:dy:ymax);    % grid point values
%
     slope = feval(filename,xm,ym);   % slope values in npts+1 by npts+1 array
%                                      %   (delta y)
     deltax = ones(npts+1);           % array of unity (delta x) values
%
%  plot slope/derivative vectors
%
     f = scale./sqrt((slope/dy).^2 + (1/dx)^2);  % scale factor for
%                                                 % constant grid length
     f2 = f/2;     % for centering vector on mesh point
%    
     xm1 = xm - deltax.*f2;
     ym1 = ym - slope.*f2;       % starting end of slope vectors
%
     xm2 = xm + deltax.*f2;
     ym2 = ym + slope.*f2;       % final end of slope vectors
%
%  There is probably a faster way of doing the following
%
     for i=1:npts+1
        for j=1:npts+1
           xp=[xm1(i,j) xm2(i,j)];
           yp=[ym1(i,j) ym2(i,j)];
	   plot(xp,yp)
        end
     end
