Matlab - interpolacja newtona

X = [-1 0 1 2];
Y = [-4 -1 0 5];
n=length(X)-1;
for i = 1:n

    w(i) = (Y(i+1) - Y(i))/(X(i+1) - X(i));
    m = length(w)-1;
    for j = 1:m
        z(j) = (w(j+1) - w(j)) / (X(j+2) - X(j));
        o = length(z)-1;
        for k = 1:o
            a(k) = (z(o+1) - z(o))/(X(o+3) - X(o));
        end
    end

end

disp(w)
disp(z)
disp(a)
x=3;
wiel = Y(1) + w(1).*(x-X(1)) + z(1).*(x-X(1)).*(x-X(2)) + a(1).*(x-X(1)).*(x-X(2)).*(x-X(3))