x=[1/2,3/4,1]
f(x)=x-cos(x)
y=f.(x)
function l(k,t)
    r=1
    for j=1:3
        if j!=k
            r*=(t-y[j])/(y[k]-y[j])
        end
    end
    return r
end
function p(t)
    s=0
    for k=1:3
        s+=x[k]*l(k,t)
    end
    return s
end
alpha=p(0)
println("alpha = ",alpha)
println("f(alpha) = ",f(alpha))
x=[x[2], x[3], alpha]
y=f.(x)
alpha=p(0)
println("alpha = ",alpha)
println("f(alpha) = ",f(alpha))
x=[x[2], x[3], alpha]
y=f.(x)
alpha=p(0)
println("alpha = ",alpha)
println("f(alpha) = ",f(alpha))

println("\nResidual errors using big numbers:")
x=big.([1/2,3/4,1])
y=f.(x)
for k=1:6
    global alpha,x,y
    alpha=p(0)
    display(f(alpha))
    x=[x[2], x[3], alpha]
    y=f.(x)
end