function pp(fs,gs,xx,yy,xlab,ylab,ttl,col,clrscr);
% function pp(fs,gs,xx,yy,xlab,ylab,ttl,col,clrscr);
%
% pp: plots phase planes for differential equations of the form
%	dx/dt=f(x,y), dy/dt=g(x,y).
%
% Usage example: 
%                pp(@f,@g,-10:10,-10:10,'x','y','ttl')
% Parameters:
%	fs, gs: user defined function
%	xx: a row vector of x-coordinates
%	yy: a row vector of y-coordinates
%	xlab,ylab,ttl: label for x-axis and y-axis, title of plot
%	col: optional string indicating the color of the vectors
%	clrscr: optional flag indicating whether df is supposed to clear
%		the screen and set 'ishold' to 1 before plotting its direction field.
%		Default is 1.

	if nargin<7
		col='b';
	end;
	if nargin<8
		clrscr=1;
	end

	[x,y]=meshgrid(xx,yy);		% create a grid of x and y-coordinates
	x=x(:)';				% then turn it into a vector
	y=y(:)';				% for plotting purposes

	xdot=feval(fs,x,y);			% compute "velocity"
	ydot=feval(gs,x,y);			% components

                                % find indices of finite (i.e., defined) slopes
	idx=find(finite(xdot) & finite(ydot));	
	xdot=xdot(idx);		% throw away undefined velocities
	ydot=ydot(idx);		% 
	x=x(idx);			% as well as corresponding 
	y=y(idx);			% coordinates

	if clrscr
		fixaxes([min(xx),max(xx),min(yy),max(yy)]);
	end

	if nargin>3
		xlabel(xlab);
	end
	if nargin>4
		ylabel(ylab);
	end
	if nargin>5
		title(ttl);
	end
    
    scale = 1;   % change to scale the plot for aesthetic purposes
    quiver(x,y,xdot,ydot,scale,col)
