function bubbles(F)
    % bubbles. Color portrait of complex-valued function F(z),
    % ex. bubbles(@sin)
    %     bubbles(@(z) .5./(z.^5-z/5) )
    
    if nargin < 1
       F = @(z)z;
    end
    
    axis(1.5*[-1 1 -1 1])
    axis square
    box on
    cla
    
    m = 256;
    colormap = [parula(m);flipud(parula(m))]; 
        
    n = 25;
    s = -1:2/(n-1):1;
    [x,y] = meshgrid(s);
    z = x + y*1i;
    circle = exp((0:32)/16*pi*1i)/n;
    
    w = F(z);
    r = abs(w);
    theta = angle(w)+pi;
    
    scale = 20;
    for k = 1:n
        for j = 1:n
            t = m*theta(k,j)/pi/scale;
            idx = ceil(scale*t+realmin);
            color = colormap(idx,:);
            p = w(k,j) + r(k,j)*circle;
            patch(real(p),imag(p),color)
        end
    end
 
    titleF = char(F);
    if (titleF(1) == '@')
        titleF(1:4) = [];
    end
    title(titleF)
    snapnow
end