function draw_dirfield(f, range, ics)
%Usage: draw_dirfield(f, range, ics)
%Simple command to draw a direction field for an ODE u'=f(t,u).
%Here "range" is a vector of the form [low, highx, lowy, highy] to
%delineate the range of the plot. The optional n x 2 array "ics" contains
%initial conditions through which solution should be drawn, forward and
%backward in time.

N = 20; %Grid size for mesh
dx = (range(2)-range(1))/N; %horizontal spacing
dy = (range(4)-range(3))/N; %vertical spacing

%Matlab's built in meshgrid to set up grid.
[t,u] = meshgrid(range(1):dx:range(2),range(3):dy:range(4));

%Compute slope at each grid point, then normalize that all vectors will
%have the same length.
s = f(t,u); %Vector <1,f(t,u)>
sn = 1./sqrt(s.^2+1); %Rescale to <sn, sn*f(t,u)>, where sn = 1/sqrt(1+f^2(t,u))
s = s.*sn;

%Direction field plot, scaled so arrow length about 80 percent of gridsize
scale = 1;
quiver(t,u,sn,s,0.8*scale); axis tight

%If applicable, sketch solutions through each initial condition, forward
%and backward in time.
if exist('ics','var')
    [n,~] = size(ics);
    hold on
    R = range(2)-range(1);
    for i=1:n
        t0 = ics(i,1);
        u0 = ics(i,2);
        if abs(t0-range(2)) > (1.0e-3)*R
            [tm,um] = ode45(f,[t0,range(2)], u0); %Forward in time
            plot(tm,um,'-k')
        end
        if abs(t0-range(1)) > (1.0e-3)*R
            [tm,um] = ode45(f,[t0,range(1)], u0); %Backward in time
            plot(tm,um,'-k')
        end
    end
    hold off
    axis(range)
end