function ans=dirfield(deriv,ax,bx,ay,by,ns)
%
% function dirfield(deriv,ax,bx,ay,by,ns)
%
% Construct the direction field for the differential equation
%        y' = deriv(x,y)
% The direction field is to be constructed on the rectangle
% [ax,bx]X[ay,by].  You are to have a program named 'deriv.m'
% with a first statement 
%        function ans=deriv(x,y)
% It should allow for the possibility that x and y are arrays
% of the same order.
% The parameter ns gives the number of subdivisions to be used
% for [ax,bx] and [ay,by]. Use ns=10 to begin. 
%
%  Reference:
%
%    Kendall Atkinson, Weimin Han,
%    Elementary Numerical Analysis,
%    Wiley, 2004,
%    ISBN: 0471433373,
%    LC: QA297.A83.2004.
%
axis([ax,bx,ay,by]);
hold

% Create the mesh at which to construct the direction vectors.
hx=(bx-ax)/ns; 
hy=(by-ay)/ns;
delta=hx/5;
x=ax+hx*[1:ns]-0.5*hx;
y=ay+hy*[1:ns]-0.5*hy;

% Use the vectors x and y to construct matrices xx and yy which
% are the coordinates of the cross-product set for x and y.
[xx, yy]=meshgrid(x,y);

% Evaluate the differential equation at (xx,yy).
d=feval(deriv,xx,yy);

% Produce the direction field using the program 'linekea'.
for i=1:ns
  for j=1:ns
    linekea(xx(i,j),yy(i,j),d(i,j),delta);
  end
end

% Adjust the scales to have them be equal for the x and y variables.
axis('equal')
hold off

function ans=linekea(x,y,dy,h)
hx=h/sqrt(1+dy^2);
hy=dy*hx;
xl=x-.5*hx; xr=x+.5*hx;
yl=y-.5*hy; yr=y+.5*hy;
plot([xl xr],[yl yr],'b')
