% Dado z(n), calcula z(n+1)
function y = cIterate(x)
  if(x == 0)
    y = 0;
  else 
    y = x - (x^4-1) / (4*x^3);
  end;
 
endfunction;

% Dado un z inicial y un valor n natural, itera el algoritmo n veces.
function y = cSolve(x, n)
  y = x;
  for i=1:n;
    y = cIterate(y);
  end;
endfunction;

% Grafica el fractal, haciendo una grilla de (2 * s + 1) x (2 * s + 1) valores
function y = fractal(s);
  y = 0;
  A = zeros(2 * s + 1, 2 * s + 1);
  for i = -s:s;
    for j = -s:s;
      A(i + s + 1, j + s + 1) = color((i + j * %i) / s);
    end;
  end;
  Matplot(A);
endfunction;

% Para un z complejo, devuelve de que color debe ir pintado, dependiendo de la proximidad a las raíces del polinomio después de 10 iteraciones.
% Si z(10) no está cerca de ninguna de las raíces, devuelve 1.
function y = color(z)
  y = 1;
  z = cSolve(z, 10);
  if(abs(z - %i) < 0.01)
    y = 2;
  end;
  if(abs(z + %i) < 0.01)
    y = 3;
  end;
  if(abs(z - 1) < 0.01)
    y = 4;
  end;
  if(abs(z + 1) < 0.01)
    y = 5;
  end;
endfunction;
  
