function phaseplane(f, xmin, xmax, ymin, ymax)
%The function phaseplane draws the direction field for the autonomous 
%system of first-order differential equations [dx/dt, dy/dt] = f(x, y). 
%The function f must be defined in a separate function file. 
%phaseplane will also generate numerical solutions to the system of 
%differential equation [dx/dt, dy/dt] = f(x, y) passing through 
%user-defined data points.
%Inputs:
%The name of the file defining the right-hand side of the system of 
%differential equations.
%Minimum value of variable x
%Maximum value of variable x
%Minimum value of variable y
%Maximum value of variable y
delta_x = (xmax - xmin)/20;
delta_y = (ymax - ymin)/20;
x = xmin:delta_x:xmax;
y = ymin:delta_y:ymax;
[xm, ym] = meshgrid(x, y);
u = zeros(size(xm));
v = u;
for i=1:length(x)
    for j=1:length(y)
        xpyp = feval(f, 0, [xm(i,j), ym(i,j)]);
        u(i, j) = xpyp(1);
        v(i, j) = xpyp(2);
    end
end
quiver(xm, ym, u, v, 'k')
axis([xmin xmax ymin ymax])
xlabel('x')
ylabel('x''')
hold on
xy0 = input('Enter initial data point [x0,y0]. Just hit the Enter key with no input to stop the program. ');
while length(xy0)~=0
    [t, z] = ode45(f, [0, 100], xy0);
    plot(z(:,1),z(:,2),'-b')
    axis([xmin xmax ymin ymax])
    xy0 = input('Enter initial data point [x0,y0]. Just hit the Enter key with no input to stop the program. ');
end