# Problem with multiplying functions to create Lagrange interpolating polynomial basis

I am trying to write some functions to generate a Lagrange interpolating polynomial basis. That is, I want to create a vector where each element is a lagrange polynomial basis function. I have some code below but it does not work. And I am not sure how to set this up to work correctly. I know that there is a package in julia for Lagrange polynomials, but I wanted to write this myself so that I understand them better.

what I would like is to have some vector of functions such that each element corresponds to a lagrange basis function. So if the vector is \psi, then \psi1should apply the first lagrange basis function at the value ofx. Here is some code I am working with, but multiplying the functions to create the polynomial basis functions is proving harder than I through. Some subtle errors in how to handle functions and multiply them, etc. I was hoping someone could help. Here is some code that I am starting with.

function lagrange_polynomial(x, i, points)
p(x) = 1;
N = length(points)
for k=1:N
if k != i
p = x -> p(x)*(x - points[k])/(points[i] - points[k])
@show p(x)
end
end
return p
end

function lagrange_polynomials_01(x, N)
h = 1.0/N

points = [(i-1)*h for i=1:N+1]
ψ_array = Array{Function}(undef, N+1);
for i=1:N+1
ψ_array[i] = x -> lagrange_polynomial(x, i, points)
end
return ψ_array, points
end


So I should be able to do something like.

julia>psi, points = lagrange_polynomials_01(0.4, 2)
julia>psi(0.4)


But then I get the error:

ERROR: StackOverflowError:
Stacktrace:
 (::var"#9#11"{Int64, Vector{Float64}, Int64})(x::Float64) (repeats 10823 times)
@ Main /media/hayagriva/Dropbox/sandbox/julia/julia_pde/julia_finite_elements/langetangen/approx1d.jl:145
 macro expansion
@ ./show.jl:1128 [inlined]
 lagrange_polynomial(x::Float64, i::Int64, points::Vector{Float64})
@ Main /media/hayagriva/Dropbox/sandbox/julia/julia_pde/julia_finite_elements/langetangen/approx1d.jl:146


So not sure how to make this run smoothly. Seems like there are some simple function tricks that I am missing. Any help is appreciated.

julia> lagrange_polynomial(x, i, points) = prod((x-points[k])/(points[i] - points[k]) for k in eachindex(points) if k ≠ i)
lagrange_polynomial (generic function with 1 method)

julia> points = 0:10
0:10

julia> lagrange_polynomials = [x -> lagrange_polynomial(x, i, points) for i in eachindex(points)]
11-element Vector{var"#42#44"{Int64}}:
#42 (generic function with 1 method)
#42 (generic function with 1 method)
⋮
#42 (generic function with 1 method)

julia> lagrange_polynomials.(0:10)
11-element Vector{Float64}:
-0.0
0.0
1.0
-0.0
0.0
-0.0
0.0
-0.0
0.0
-0.0
0.0
`