function draw_auto_phase(f1, f2, viewrange, ics)

xh = viewrange(2)-viewrange(1);
yh = viewrange(4)-viewrange(3);
N = 20;
hx = xh/N; hy = yh/N;
[X1,X2] = meshgrid(viewrange(1):hx:viewrange(2),viewrange(3):hy:viewrange(4));

v1 = f1(X1,X2);
v2 = f2(X1,X2);
len = sqrt(v1.^2 + v2.^2); %Rescale arrows to all have same length---typically looks better
v1 = v1./len;
v2 = v2./len;

quiver(X1,X2,v1,v2,0.5,'-b'); %0.5 defines arrow lengths---adjust to taste.
axis(viewrange)

syms p1 p2;
sol = solve([f1(p1,p2)==0,f2(p1,p2)==0],[p1,p2]);
n = length(sol.p1);
hold on
scatter(sol.p1,sol.p2,'filled','MarkerFaceColor',[0 0 0]);

%If applicable, sketch solutions through each initial condition, forward
%and backward in time.
if exist('ics','var')
    f = @(t,x) [f1(x(1),x(2)) f2(x(1),x(2))]';
    [n,~] = size(ics);
    R = 5; %Will solve for -R < t < R
    for i=1:n
        %Assume ics are a n by 2 matrix
        u0 = ics(i,:);
        [~,um] = ode45(f,[0,R], u0,odeset('RelTol',1e-4)); %Forward in time
        plot(um(:,1)',um(:,2)','-r')

        [~,um] = ode45(f,[0,-R], u0,odeset('RelTol',1e-4)); %Backward in time
        plot(um(:,1)',um(:,2)','-r')
    end
    hold off
end
hold off
drawnow