clear
epsilon = 0.01;
a = [1 4 -72 -214 1127 1602 -5040];
color = [0 0 0; 255 0 0; 255 255 0; 0 255 0; 0 0 255; 255 0 255]/255;
grad = [fliplr(0:0.1:1) 0:0.1:1];
xlim = [-9 8];
ylim = [-2000 2000];
x0 = 10;
x = [xlim(1):.01:xlim(2)];
roots(1) = newton(a,x0,epsilon);
b = a;
for i = 2:length(a)-1
[y a] = horner(b(i-1,:),roots(i-1));
b(i,:) = [0 a];
roots(i) = newton(b(i,:),roots(i-1),epsilon);
endfor
b(length(a),:) = b(1,:);
for i = 1:length(a)
# fancy graphics
for j = 1:length(grad)
shade = grad(j)*([1 1 1]-color(i,:));
hold off
plot(x,polyval(b(i,:),x),'color',color(i,:)+shade,'linewidth',3)
hold on
plot(x,polyval(b(1,:),x),'color',color(1,:),'linewidth',3)
plot(x,zeros(size(x)),'--k','linewidth',3)
for k = 1:i-1
plot(roots(k),0,'o','color',color(k,:),'markersize',1,'linewidth',3)
endfor
if j < length(grad)/2
plot(roots(i),0,'o','color',color(i,:)+shade,'markersize',1,'linewidth',3)
else
plot(roots(i),0,'o','color',color(i,:),'markersize',1,'linewidth',3)
endif
axis([xlim ylim])
print(strcat("frame",num2str(j+length(grad)*(i-1)),".eps"))
endfor
endfor

function z = newton(a,x0,epsilon)
x1 = epsilon*2+x0;
loops = 0;
for i = 1:length(a)-1
b(i) = a(i)*(length(a)-i);
endfor
while abs(x0-x1) > epsilon && loops < 500
x0 = x1;
f = horner(a,x0);
fp = horner(b,x0);
x1 = x0 - f/fp;
loops++;
endwhile
z = x1;
endfunction

function [y b] = horner(a,x)
b(1) = a(1);
for i = 2:length(a)
b(i) = a(i)+x*b(i-1);
endfor
y = b(length(a));
b = b(1:length(b)-1);
endfunction