a = [-1,1,-1,1,-1,1,1]; pow = [1:length(a)]-1;
b = a;  c = a; % initialize b & c, just good programming practice
p = @(x) sum((x.^pow).*a); 
dp = @(x) sum((x.^(pow(2:end)-1)).*a(2:end).*pow(2:end));

x0 = .5;
i = length(a);
while i > 2
    i = i-1;
    b(i) = a(i) + b(i+1)*x0;
    c(i) = b(i) + c(i+1)*x0;
end
b(1) = a(1) + b(2)*x0;
c(2)-dp(x0)
b(1)-p(x0)

x = linspace(-1,1,100); y = x;
for i = 1:100
    y(i)=p(x(i));
end
hold off; plot(x,y)
title('A plot of P(x)')

p0 = .1;
for j = 1:20;
    i = length(a);
    while i > 1
        i = i-1;
        b(i) = a(i) + b(i+1)*p0;
        c(i) = b(i) + c(i+1)*p0;
    end
    p0 = p0 - b(1)/c(2);
end
p0

plot(x,y); hold on;
title('A plot of P(x)')
plot(x,0*x) % plot the x-axis
plot(p0,0,'ro') % plot the zero to see that we are computing it