function nbasin(func, deriv, lims, res, iter)
%NBASIN Graphs the basins of attraction for Newton's method
%   NBASIN(FUNC, DERIV, LIMS, RES, ITER) performs Newton's method
%   for a grid of initial guesses in the complex plane and
%   color-codes them according to where they end up after ITER
%   iterations.  Formulas for both the function FUNC to find the
%   roots of and its derivative DERIV must be given as strings.
%   The range of real and imaginary parts must be given as a
%   vector LIMS of 4 numbers, and RES is the resolution of the
%   graph (the number of grid points in each direction).
%
%   Example:
%   nbasin('x^3 - x + 0.7', '3*x^2 - 1', [-5 5 -5 5], 300, 50)
%
%   The algorithm for the color-coding is crude and chooses decent
%   colors for the example above, but will probably choose colors
%   poorly for other functions.  After ITER iterations from a
%   given initial point, the real and imaginary parts X and Y of
%   the result determine the red, green, and blue intensities at
%   that point according to three formulas in the M-file nbasin.m.
%   Changing these formulas may improve the coloring.  If the
%   contrast between colors is too low, increasing ITER may help.

f = inline(vectorize(func));
df = inline(vectorize(deriv));
x = linspace(lims(1), lims(2), res);
y = linspace(lims(3), lims(4), res);
[X, Y] = meshgrid(x, y);
Z = X + i*Y;
for k = 1:iter
    Z = Z - f(Z)./df(Z);
end
X = real(Z);
Y = imag(Z);

red = X;
green = Y;
blue = -Y.^2;

mat = zeros(res, res, 3);
rmin = min(min(red));
rmax = max(max(red));
mat(:,:,1) = (red - rmin)/(rmax - rmin);
gmin = min(min(green));
gmax = max(max(green));
mat(:,:,2) = (green - gmin)/(gmax - gmin);
bmin = min(min(blue));
bmax = max(max(blue));
mat(:,:,3) = (blue - bmin)/(bmax - bmin);
image(x, y, mat)
set(gca, 'YDir', 'normal')
axis equal tight
