function julijaiter(rez, iter,C)
% julijaiter(10000, 500,0.5i)

pts = zeros(rez,2);
x = rand(1);
y = rand(1);
a=real(C);
b=imag(C);

f=@(x,y) sqrt((x-a+sqrt((x-a)^2+(y-b)^2))/2);
for j = 1:rez+iter
    p1 = rand(1);    
        if p1 < 3/4        
            xn = f(x,y);
            yn = (y-b)/(2*f(x,y));
        else
            xn = -f(x,y);
            yn = -(y-b)/(2*f(x,y));
        end        
        x = xn;
        y = yn;
    pts(j,1) = x;
    pts(j,2) = y;
end%j

xs = pts(iter:end,1);
ys = pts(iter:end,2);
plot(xs,ys,'.g')
axis equal