function M=julia(zmin,zmax,hpx,k,c)
% c e' il parametro
% k il numero di iterazioni
% hpx e' il numero di pixel orizzontali

  vpx=round(hpx*abs(imag(zmax-zmin)/real(zmax-zmin)));
  [zRe,zIm]=meshgrid(linspace(real(zmin),real(zmax),hpx),
                     linspace(imag(zmin),imag(zmax),vpx));
  z=zRe+i*zIm;
  M=zeros(vpx,hpx);
  
  for s=1:k
    m=abs(z)<2;
    M(m)=M(m)+1;
    z(m)=z(m).^2+c;
  endfor
  M(m)=0
endfunction
