function richtungsfeld(df, xlim, ylim, n )
if nargin < 4, n = 40; end

gray = [127, 127, 127] / 255;

xpts = linspace(xlim(1), xlim(2), n);
ypts = linspace(ylim(1), ylim(2), n);

[x, y] = meshgrid(xpts, ypts);

px = ones(size(x)); py = df(x, y);
nm = sqrt(px.^2 + py.^2);
px = px ./ nm; py = py ./ nm;

quiver(x, y, px, py, 0.75, 'LineWidth', 1.0, 'Color', gray);
axis([xlim ylim])
end