using Plots
runge = x->1/(1+x^2)
f = x->1/(1+25x^2)
#P2D to macierz Nx2 o wierszach [x y]
# współczynniki ilorazów dla interp. Newtona
function _diff(P2D)
    b, n = Array(P2D[:, 2]), length(P2D[:, 1])
    for  i = 2:n
        for j = n:-1:i
            b[j] = (b[j] - b[j-1])/(P2D[j, 1] - P2D[j-i+1, 1])
        end
    end
    return b
end

# Zwraca wartość ilorazu.
function diff(P2D)
    last(_diff(P2D))
end

function interp_newton(P2D)
    return function(x)
        b, n = _diff(P2D), length(P2D[:, 1])-1
        sum, p = b[1], 1
        for i in 1:n
            p *= (x - P2D[i, 1])
            sum += b[i+1]*p
            end # nie chce mi się robić Hornera
        sum
    end
end;

err = f_interp->x->abs(f(x)-f_interp(x))
w9_eq = interp_newton(reduce(vcat, map(x->[x f(x)], collect(range(-1,stop=1, length=10)))))             # w równoodległych
w9_T10_zeros = interp_newton(reduce(vcat, map(x->[x f(x)], map(j->cos(pi*((2*j-1)/(2*10))), [1:10;])))) # w zerach Czebyszewa
equal_plot = begin
    plot(f, -1, 1, line=(color=:blue), label="f", title="equal")
    plot!(w9_eq, line=(color=:orange), label="equal")
    plot!(err(w9_eq), fill=(0, 0.1, :red), label="EQ error")
end
Chebyshev_zeros_plot = begin
    plot(f, -1, 1, line=(color=:blue), label="f", title="Chebyshev's zeros")
    plot!(w9_T10_zeros, line=(color=:orange), label="Ch. zeros")
    plot!(err(w9_T10_zeros), fill=(0, 0.1, :red), label="CH. zeros error")
end
row_plot = plot(equal_plot, Chebyshev_zeros_plot, layout=(1,2), fmt=:png, size=(900, 300), ylims=(-0.5, 2))