function [] = state_dfield(xylimits,npts,filename,scale)
%  Like dfield.m, except plots derivative values in state or phase space.
%
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
mesh = [xm,ym];
mesh = reshape(mesh,npts+1,npts+1,2); %convert to  array (2 state variables)
t = 1;    %  dummy variable for getting derivative
     for i=1:npts+1
        for j=1:npts+1
           z=mesh(i,j,:);       %  state variables for the mesh point
           zdot = feval(filename,t,z);  %  derivatives of the state variables
           z1dot(i,j) = zdot(1);
           z2dot(i,j) = zdot(2);
        end
     end
deltax = ones(npts+1);  %  array of unity (delta x) values
slope = z2dot./z1dot;   %  slope values in 2D space
%
%  Plot vectors with constant length
%
f = scale./sqrt((slope/dy).^2 + (1/dx)^2);  %  scale factor
f2 = f/2;  %  for centering on mesh point
%
xm1 = xm - deltax.*f2;   %  start of slope vectors
ym1 = ym - slope.*f2;
xm2 = xm + deltax.*f2;
ym2 = ym + slope.*f2;    %  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


