function [Z,renderMat]=fractal(n,d,l)
%% Z is the n-by-n output of solutions to the equation ((x+i*y)^d)-1=0
%% n is the number of points in the domain.
%% d is the degree of the polynomial
%% l is the length of the real and imaginary axes used.
	clc
	maxIter=20;					% Maximum number of Newton iterations to be done
	cx=0;						% Center of domain (Real axis)
	cy=0;						% Center of domain (Imag axis)
	tol = 1e-10;				% Tolerance to determine closeness to 0.0
	x=linspace(cx-l,cx+l,n);	% Real space
	y=linspace(cy-l,cy+l,n);	% Imaginary space
	[X,Y]=meshgrid(x,y);
	Z=X+i*Y;
	
	%% Perform Newton iterations
	for k=1:maxIter;
		Z=Z-(f(Z,d)./fprime(Z,d));
	end

	%%%%%%%%%%%%%%%%%%%%%%%%%%
	%% Plotting the fractal %%
	renderMat=zeros(n,n);			% Matrix to render each root convergence as its own color

	%% Find d roots of unity, and the jth mask
	for j=1:d
		root=exp(2*pi*i/d)^j;		% The jth root of the d roots of unity
		Mj=abs(Z-root);				% Where did Z converge close to root
		mask=(Mj<=tol)*j;			% Each root gets a unique number in [1,d]
		renderMat=renderMat+mask;	% Add it to the rendering matrix
	end
	colormap(hsv);		% Set the color map
	imagesc(renderMat)	% Render the fractal
	shading flat;		% Turn off grid lines
	axis('square','equal','off');
	imgStr=sprintf('newton_%i_d_%i_L_%i.png', n,d,l);
	titleStr=sprintf('Newton Method: N=%i Tol:%3.2e Deg:%i Len:%i', n,tol,d,l);
	title(titleStr)
	print('-dpng',imgStr)	% Output image to file

end

function y=f(x,d);
	y=(x.^d)-1;
end

function y=fprime(x,d);
	y=d*(x.^(d-1));
end
