function [zz] = fun_brot(x_min,x_max,y_min,y_max)
%

% set parameters and initial conditions
xy_res = 1000;
k_max = 500;
thresh = 2.0*2.0;
% create grid
[xx, yy] = meshgrid([x_min:(x_max-x_min)/xy_res:x_max],[y_min:(y_max-y_min)/xy_res:y_max]);
% turn into vector format
n_max = numel(xx);
v(:,1) = reshape(xx,[n_max,1]);
v(:,2) = reshape(yy,[n_max,1]);
v(:,3) = zeros(n_max,1);

% solve equation
for n=1:n_max
    k = 0;
    x = v(n,1);
    y = v(n,2);
    a = 0.0;
    b = 0.0;
    while ( (k <= k_max) && ((a*a + b*b) <= thresh) ),
        k = k+1;
        a_tmp = a;        b_tmp = b;
        a = a_tmp*a_tmp - b_tmp*b_tmp + x;
        b = 2*a_tmp*b_tmp + y;
    end
    if (k <= k_max),
       v(n,3) = (1.0 - k/k_max)^3;
    else
       v(n,3) = 0.0;
    end
end

% return results as a matrix
zz = reshape(v(:,3),[xy_res+1,xy_res+1]);
