function vectorfield(f,y1val,y2val,t)
% vectorfield(f, a1:dy1:b1, a2:dy2:b2)
%   
%   plot direction field for first order ODE y' = f(t,y)
%   using y1-values from a1 to b1 with spacing of dy1
%   using y2-values from a2 to b2 with spacing of dy2
%
%   f is defined as f = @(t,y) [ ... ; ... ]
%     where "..." are expressions with t, y(1), y(2)
%
%   Usually this is used for an autonomous ODE.
%   If f depends on t the vector field is plotted for t=0
%   (a t-value can be specified as an optional last argument).
%
% Example: y1' = y2
%          y2' = -sin(y1)
%   Show vector field for y1 in [-2,8], y2 in [-2.5,2.5], use
%   spacing of .5 for y1, .25 for y2
%
%   f = @(t,y) [ y(2) ; -sin(y(1)) ]; 
%   vectorfield(f, -2:.5:8, -2.5:.25:2.5)

if nargin==3
  t = 0;
end
col = [1 .5 .5];              % color for vector field
[y1m,y2m] = meshgrid(y1val,y2val);
dy1 = y1val(2) - y1val(1); 
dy2 = y2val(2) - y2val(1);
yp1 = arrayfun(@(y1,y2)[1 0]*f(t,[y1;y2]),y1m,y2m); 
yp2 = arrayfun(@(y1,y2)[0 1]*f(t,[y1;y2]),y1m,y2m); 
quiver(y1val,y2val,yp1,yp2,'Color',col); 
axis([y1val(1)-dy1/2,y1val(end)+dy1/2,y2val(1)-dy2/2,y2val(end)+dy2/2])

