function [J,x,y]=juliaset(area,Nmax,N,a,b,c)

x=linspace(area(1),area(2),N);
y=linspace(area(3),area(4),N);
[X,Y]=meshgrid(x,y);
Z=X+1i*Y;

limit=grense(a,b,c);
Z0=Z/sqrt((area(2)^2+area(4)^2));
J=zeros(N,N);
for i=1:Nmax,
    K=find(abs(Z)<=limit);
    Z(K)=poly(Z(K),a,b,c);
    J(K)=J(K)+1;
end;
J=(J/Nmax+1./(1+abs(Z)));

%-------------------
function p=poly(z,a,b,c)
p=a*z.^2+b*z+c;

%-------------------
function l=grense(a,b,c)
aa=abs(a); ab=abs(b); ac=abs(c);
l=0.5*(ab+1+sqrt((ab+1)^2+4*aa*ac))/(aa);