From: https://thoughtforms.life/halleys-method-fractal-art/

newtonsmethod.m

% Copied from thoughtforms.life
function Newtonsmethod
global cancelled;
lambda = 1;
cancelled = 0;
r=0;
e = 2.71828182845904523;
epsilon = 0.00001;
% make the default colormap.
clrmp = [0 0 0; 0 0 1; 0 1 0; 1 0 0; 0 1 1; 1 0 1; 1 1 0; 1 1 1];
fname = input('File to save movie in: ','s');
steps = 0;
r_begin = 0;
r_end = 0;
% (sum(r_end)==0) || sum(r_begin)==0 || 
while (r_end<=r_begin)
    r_begin = input('Enter begin value for real axis: ');
    r_end = input ('Enter end value for real axis: ');
end
i_begin = 0;
i_end = 0;
% (sum(i_end)==0) || sum(i_begin)==0 || 
while (i_end<=i_begin)
    i_begin = input('Enter begin value for imaginary axis: ');
    i_end = input ('Enter end value for imaginary axis: ');
end
while (steps<32)
    steps = input('Enter resolution in pixels: ');
end
i_step = (i_end-i_begin)/steps;
r_step = (r_end-r_begin)/steps;
formula=input('Enter formula involving the complex variable z and movie parameter t: ','s');
derivative=diff(formula,'z');
q1 = sprintf('First derivative is %s',char(derivative));
disp(q1);
% sanity check for 1st derivative
% t=1.0; 
% z = 10.0; 
% zy1 = eval(formula);
% z = 10.0+epsilon;
% zy2 = eval(formula);
% slope = (zy2-zy1)/epsilon;
% z = 10.0+epsilon/2.0;
% zz = eval(derivative);
% if (abs(zz-slope)>epsilon)
%     disp('That does not seem like the right first derivative!'); 
% end
derivative2 = diff(derivative,'z');
q2 = sprintf('Second derivative is %s',char(derivative2));
disp (q2);
% sanity check for 2nd derivative
% t=1.0; 
% z = 10.0; 
% zy1 = eval(derivative);
% z = 10.0+epsilon;
% zy2 = eval(derivative);
% slope = (zy2-zy1)/epsilon;
% z = 10.0+epsilon/2.0;
% zz = eval(derivative2);
% if (abs(zz-slope)>epsilon)
%     disp('That does not seem like the right second derivative!'); 
% end
movie_begin = 0;
movie_end = 0;
movie_step = input('How many frames for movie? ');
if (sum(movie_step)==0 || movie_step<2)
    if (length(fname)>0)
        filename = sprintf('%s.jpg',fname);
        message = sprintf('Will make single frame into file %s.', filename);
        disp(message);
    else
        disp('Will make single frame. ');
    end
    movie_begin_=0;
    movie_end=1;
    movie_step=1;
    message = sprintf('Estimated completion time = %g hours.', ...
              0.00481004 * (steps*steps*movie_step) /(60*60)); 
    disp(message);
else
    if (length(fname)>0)
        filename = sprintf('%s.avi',fname);
        message = sprintf('Will make movie into file %s, estimated size = %d Kbytes.', ...
                  filename, 0.984*movie_step*steps*steps/1000);
        disp(message);
    else
        disp('Will make movie. ');
    end
    message = sprintf('Estimated completion time = %g hours.', ...
              0.00481004 * (steps*steps*movie_step) /(60*60)); 
    disp(message);
    while (movie_end<=movie_begin)
        movie_begin = input('Enter begin value for t: ');
        movie_end = input ('Enter end value for t: ');
    end
end
if (movie_step<0)
    movie_step=1;
end
clrs = 0;
while (clrs<1)
   clrs = input('Enter colormap to use: 1 = old style (8), N = smooth: '); 
end
if (length(fname)>0)
    datafile = sprintf('%s.txt',fname);
    [fid message] = fopen(datafile,'w');
    if (fid==(-1))
        disp(message)
    else
        fprintf(fid,'minx=%g, maxx=%g, miny=%g, maxy=%g, \nt is from %g to %g, \nformua = %s.\n', ...
            r_begin,r_end,i_begin,i_end,movie_begin,movie_end,formula);
        fprintf(fid,'resolution = %d, frames = %d, method = Halley. \n',steps,movie_step);
    end
else
    fid = (-1);
end
if (clrs==1)
    colormap(clrmp);
    clrs = 8;
else
    colormap(hsv(clrs));
end
cmap = colormap;
if (movie_step>1 && length(fname)>0)
    aviobj=avifile(filename,'colormap',cmap,'compression','None','fps',7);
end
if (clrs==8)
    klimit=15; % old-style algorithm
else
    klimit=3*clrs; % for bigger colormap; probably takes longer!
end
begint = clock;
w_movie = waitbar(0,'Movie being made...','CreateCancelBtn',@cancel_callback);
% outer movie loop - each one is a frame of the movie
frames_done = 0;
B = zeros(steps+1,steps+1, 'uint8'); % static init for performance gain
axx = r_begin : r_step : r_end;
axy = i_begin: i_step: i_end;
stats = zeros(1,51,'uint16');
for t = movie_begin : ((movie_end - movie_begin)/movie_step) : movie_end
    frames_done = frames_done + 1;
    time1 = now;
    % w_real = waitbar(0,'Real component...');
    rr = 1;
    for (r=r_begin : r_step : r_end)
        if (cancelled==1 || not(ishandle(w_movie)))
            break;
        end
        % plot settings
        plothandle = imagesc(axx,axy,B);
        ctitle = sprintf('Halley, formula = %s, t = %g ',formula,t);
        title(ctitle);
        drawnow;
        % waitbar((r-r_begin)/(r_end-r_begin),w_real);
        % w_imag = waitbar(0,'Imaginary component...');
        ii=1;
        for (i=i_begin: i_step: i_end)
            % waitbar((i-i_begin)/(i_end-i_begin),w_imag);
            % generate image
            z = complex(r,i);
            for (k=1:50)
                z = z - eval(formula)/eval(derivative) / (1 - ...
                    (eval(formula)*eval(derivative2))/(2*power(eval(derivative),2.0)));
                g = power(real(z),2)+power(imag(z),2);
                if (k>1 && abs(g-gprev)<epsilon)
                    break;
                end
                gprev = g;
            end
            % stats(k) = stats(k) + 1;
            if (abs(g-gprev)<0.00001 && mod(k,2)==0)
                B(rr,ii) = int8(mod(k,clrs));
            else
                B(rr,ii) = 0;
            end
            ii = ii + 1;
        end
        % close(w_imag);
        rr = rr + 1;
    end
    % delete(w_real);
    if (cancelled==1 || not(ishandle(w_movie)))
        break;
    end
    % plot settings
    plothandle = imagesc(axx,axy,B);
    ctitle = sprintf('Halley, formula = %s, t = %g ',formula,t);
    title(ctitle);
    drawnow;
    if (movie_step>1 && length(fname)>0)
        % add the frame to the movie file
    	% F = getframe(gca);
        F = im2frame(B,cmap);
        aviobj = addframe(aviobj,F);
    elseif (length(fname)>0) 
        imwrite(B,colormap,filename,'jpg');
    end
    time2 = now;
    timetook = time2 - time1;
    time_remaining = timetook * (movie_step - frames_done);
    estimate_done = datestr(now + time_remaining);
    message = sprintf('Predicted endpoint = %s',estimate_done);
    waitbar((t-movie_begin)/(movie_end-movie_begin),w_movie,message);
    if (movie_step==1)
        break;
    end
end
endt = clock;
% close the files
if (movie_step>1 && length(fname)>0)
    aviobj = close(aviobj);
end
if (fid==(-1))
    disp('Warning: no data file created - remember the parameters!');
else
    e = etime(endt,begint);
    fprintf(fid,'Computation took %g minutes, or %g seconds/pixel. \n', ...
            e/60.0,e/(steps*steps*movie_step));
    % this should really convert to hours etc.
    fclose(fid);
end
disp('All done. ');
if (ishandle(w_movie)) % otherwise it got cancelled and is already closed
   delete(w_movie);
end
end

function n=reverse(z)
rr = imag(z);
ii = real(z);
n = complex(rr,ii);
end

function cancel_callback(h,varargin)
fig_handle=get(h,'parent');
delete(fig_handle);
cancelled = 1;
end