function dirfieldsys2(fx_func_handle_01,fy_func_handle_02,xspan_vector_03,yspan_vector_04)
% dirfieldsys2(fx,fy, xspan, yspan)
%   
%    Plots direction field for first order 2x2 system.
%    fx and fy are vectorized @ functions in x and y, 
%    xspan = [xmin xmax] and yspan = [ymin ymax]; 
%    give bounds for x and y.

% make sure that we get any variables we need from the workspace
    datafile = [ tempname(cd), '.mat' ];
    evalin( 'base', ['save(''', datafile, ''')'] );
    load(datafile);
    % and delete the temporary file
    delete(datafile)
    
% ensure that we're using element-by-element operations
    fxstr = char( fx_func_handle_01 );
    fxstr = strrep( strrep( fxstr, '^', '.^' ), '..^', '.^' );
    fxstr = strrep( strrep( fxstr, '*', '.*' ), '..*', '.*' );
    fxstr = strrep( strrep( fxstr, '/', './' ), '../', './' );
    % I'd rather use str2func here, but that doesn't have access
    % to any variables that we've imported from the workspace
    fxint = eval( fxstr );

    fystr = char( fy_func_handle_02 );
    fystr = strrep( strrep( fystr, '^', '.^' ), '..^', '.^' );
    fystr = strrep( strrep( fystr, '*', '.*' ), '..*', '.*' );
    fystr = strrep( strrep( fystr, '/', './' ), '../', './' );
    fyint = eval( fystr );
    
    % better variable names
    xs_dfs2 = xspan_vector_03;
    ys_dfs2 = yspan_vector_04;
    
    num=20;                       % Number of arrows in each direction.
    xh=(xs_dfs2(2)-xs_dfs2(1))/num;   % Increments for x values.
    yh=(ys_dfs2(2)-ys_dfs2(1))/num;   % Increments for y values.
    xnum=(xs_dfs2(2)-xs_dfs2(1))/xh;  
    ynum=(ys_dfs2(2)-ys_dfs2(1))/yh;

    xpts=[xs_dfs2(1):xh:xs_dfs2(2)];  %  x-coordinates for arrows
    ypts=[ys_dfs2(1):yh:ys_dfs2(2)];  %  y-coordinates for arrows
    [xm,ym]=meshgrid(xpts,ypts);  %  meshgrid of coordinates

    xp=fxint(xm,ym);                 %  x' values
    yp=fyint(xm,ym);                 %  y' values

    % Next segment scales all of the vectors in the direction field
    % so they have the same length.

    scalefactor=1/25*(xs_dfs2(2)-xs_dfs2(1));  
    r = sqrt(xp.^2+yp.^2);
    dx = (xp./r)*scalefactor;
    dy = (yp./r)*scalefactor;

    % Direction field with square axes.

    quiver(xm,ym,dx,dy,0, 'k');    
    axis([xs_dfs2(1) xs_dfs2(2) ys_dfs2(1) ys_dfs2(2)]);

    if ((xs_dfs2(2)-xs_dfs2(1)) == (ys_dfs2(2) - ys_dfs2(1)))      
        axis square
        sq=1;
    else
        sq=0;
    end

end