Let's have a look how to implement Lagrange polynomials and interpolation with Lagrange polynomials on the computer using Python. So, first let's initialize the Lagrange polynomials. The equation again is given here and we define a function that we of course call Lagrange, and it will calculate the polynomials for order N, that we inject as a parameter. It returns basically the Lagrange polynomials at a highly resolved vector in x defined between minus one and one. For the given co-location points xi. These do not necessarily have to be the GLL points and I encourage you to play around with this, but of course for our spectral elements simulation, we initiate them using the GLL points. So, this is a fairly general routine that we can use first of all to look at what Lagrange polynomials look like. The way to do this is as you see here, we first give an order, let's say N equals two, and we initialize the space here between minus one and one with 1,000 points. We get the Gauss-Lobatto Legendre points from a function that I encourage you to look at. Then basically, we just call the Lagrange routine and plot it. That's given here for order two. So, we can see three Lagrange polynomials order two. The red line, you can see starts at one minus one with the value of one and then basically goes through the central point here, and the final endpoint also with value zero as expected. Of course, we know we can increase the order, let's go up to 10 and run the program again, and you want to return so you can see the beautiful polynomials, and you can also appreciate that away from the boundary, now the extreme values become denser and denser with points, where the extreme values become denser and denser, than we've seen this, this is an expression of these GLL points. Now, how does the interpolation work with these Lagrange polynomials. Well, we've implemented this further down and we define a function in the interval between minus one and one. As an exercise, I encourage you to change that function here, we simply use sine pi times x. So, basically this is one period that fits one wavelength that fits into that interval minus one and one. Let's be a little stupid and start with an order n equals two, and I would encourage you to go through the way this is then implemented. First, we have to basically get the Lagrange polynomials for this order which is done by this loop here. So, we initialize the Lagrange polynomials on the grid. Then, this is the central and very simple algorithm here, this loop here. We then basically simply multiply the value of our original function at the co-location points. This is the fij with the appropriate Lagrange polynomials and sum everything up. This is then the vector s and then we can compare the original function or sine function with the interpolative function. So, let's do this. Now, here's the result. Here we see the original function as the black dashed line, and the interpolating function is the green line. Oops, the green line is actually flat at zero everywhere. So, that's the indication that with order N equals two, we're not actually doing a good job interpolating our sine function. So, let's increase N to, let's say four and see what happens. We see that looks already much better with, in this case then five co-location points. The error is calculated here between, this is now the error between the continuous functions in the whole domain, and it's here at 14 percent which depending on what you want, this is not a very good approximation. Then if you go higher, let's say really high order 10, you can now see, well, this is a great result without, basically don't see a difference. There's almost 10 to the minus three percent error. So, that's really very, very small and it's a nice indication how you converge with higher order, with increasing order of your Lagrange polynomials with which you interpolate. Here the sine function you get more and more accurate. Now, I encourage you as an exercise to basically take this piece of code and invent your own function and check out how you can interpolate your arbitrary function using Lagrange polynomials.