function m  = julia(zmax,c,N)
% Compute the number of iterations required to escape the orbit of each
% point z in 500x500 mesh in [-zmax zmax]^2 in iteration z_n+1 = z_n^2 +c
% parameter c is given as argument. To avoid infinite loops, a limiting
% parameter N is given: if the escape does  not happen in N iterations, it
% is assumed not to happen. 
x = linspace(-zmax,zmax,500); 
y = x'*1i;
[x,y]=meshgrid(x,y);
Z = x+y;
m = zeros(500);
for k = 1:500
    for p = 1:500
        m(k,p) = escapeVelocity(Z(k,p),c,N);
    end
end
end

function n = escapeVelocity(z0,c,N)
% Function computes the necessary number of iterations to escape from orbit
% of z0 in Prisoner set iteration
n = 0;
zn = z0^2+c;
while(n<=N&&abs(zn)<=2)
    zn=zn^2+c;
    n=n+1;
end
end