function flow(fun,xlims,scale)
Nx = 20;
[xx,yy] = meshgrid(xlims(1):(xlims(2)-xlims(1))/Nx:xlims(2), 0);
N = length(xx);
px = zeros(1,N);
py = zeros(1,N);
for k = 1:N
    px(k) = feval(fun,xx(k));
end
if nargin <= 2
    scale = 0.3;
end
quiver(xx,yy,px,py,scale)
axis equal