function julia
% julia()
%
% Berechnet Julia-Mengen
% Parameter werdem im Programm angepasst

beispiel=2;
daempfung=0;
step = 0.01;

switch beispiel
    case 1
        % Datensatz 1
        %------------------------------------------------------------------
        f=@(x) x^3-1;
        df=@(x) 2*x^2;
        z = [-0.500000000000000 + 0.866025403784439i
             -0.500000000000000 - 0.866025403784439i
             1.000000000000000 + 0.000000000000000i];
        gamma=10;
        %------------------------------------------------------------------
    case 2
        % Datensatz 2
        %------------------------------------------------------------------
        f=@(x) x^5-x^3-1;
        df=@(x) 5*x^4-3*x^2;
        z = [1.236505703391499
             0.340794866197006 + 0.785423103049449i
             0.340794866197006 - 0.785423103049449i
             -0.959047717892756 + 0.428365956254189i
             -0.959047717892756 - 0.428365956254189i];
        gamma=10;
        %------------------------------------------------------------------
end

A=cell(length(-2:step:2),length(-2:step:2));
nsteps=size(A,1);
for j=1:nsteps
    for k=1:nsteps
        A{k,j}=newtonsol(-2+step*(j-1)+1i*step*(k-1)-2i,f,df,z,daempfung,gamma);
    end
end
A=cell2mat(A);
A=flipud(A);

figure
surf(-2:step:2,-2:step:2,A,'EdgeColor','none','FaceColor','interp')
colormap hsv
view(0,90)
axis square
end




function [col,x] = newtonsol(x0,f,df,z,daempfung,gamma)
% [col,x] = newtonsol(x0,f,df,z,daempfung,gamma)
% 
% Newton-Verfahren (zur Berechnung von Julia-Mengen)
%
% Eingabe:  x0        : Startwert
%           f         : nichtlineare Funktion
%           df        : Ableitung von f
%           z         : Liste der exakten Nullstellen
%           daempfung : falls 0, dann keine Daempfung, sonst Daempfung
%                       (Standard: 0)
%           gamma     : Schaetzung fuer Lipschitzkonstante (Standard: 10)
% Ausgabe:  col       : (Farb)-Wert, der zum Startwert x0 gehoert
%           x         : berechnete Nullstelle

if nargin<5
    daempfung=0;
end
if nargin<6
    gamma=10;
end

tol=1e-4;
nstep=650;

x=x0;
count=0;
while min(abs(z-x))>tol && count<nstep
      dfx=df(x);
      w=f(x)/dfx;
      if daempfung
          beta=abs(dfx);
          alpha=abs(w);
          lambda = min([1,1/(alpha*beta*gamma)]);
      else
          lambda=1;
      end
      x = x - lambda*w;
      count=count+1;
end

r = find(abs(z-x) < tol);

col=-1;
if ~isempty(r)
    switch r
        case 1
            col = 0+count/nstep;
        otherwise
            col = r+count/nstep;
    end
end
end