So, let's see how the spectral element method in one-dimension works with a heterogeneous example. We have to stick to one dimensions here. Extending the spectral element method to do this is way more complicated and involved than in the finite difference case. But already we see some interesting aspects here in the one-dimensional case. So, let's have a look. The initialization is very similar to the homogeneous case. Of course we use here 150 elements that's Ne. The order of our simulation is capital N equals three. So, we keep the size of the elements here constant, which is why we again divide the maximum physical domain, in this case eight kilometers, we divide this by the number of elements to obtain the length of each element. Now we have a quick look at how we initialize a low velocity zone at the center of our domain. Actually, what we do is we will only change the shear modulus Mu, the density will remain constant. So, this is done in this little window here. The number of elements that spends the low velocity zone will be 25, and we will reduce the velocity by a factor of 0.4. So, that means we have to first calculate the width of the fault zone in grid points, which is given here and of course depends on the order of our scheme and we keep this general. Then we can start initializing first the velocity field here that's the vs, is a shear velocity. Of course N times Ne plus ONE, you remember that. This is the global number of co-location points that are its number of degrees of freedom for all elements of our parameters which is shear velocity, density, and also Mu. So, we initialize the velocity, reducing it. Well, actually we reduce it by a factor of 0.4. So, it's actually 60 percent slower the velocity. Then we initialize Mu here then this is now a heterogeneous model. Just note again Mu, Rho, and vs are all vectors. Now, the mass matrix will be calculated in exactly the same way. Actually here, Rho I just mentioned it is a vector so Rho could also change and that would make the mass matrix elements also vary at each co-location point but in this case Rho was constant. So, the mass matrix is given here, just a part of the mass matrix is given here, the central part again is diagonal. In fact, you see the green elements are the ones at the element boundaries because they have different values. The color scale is not important here but because they actually contain the sum of two elements of adjacent finite elements. Now, the stiffness matrix is initialized here and what is changing here as you can see that compared to the homogeneous case, now Mu has to be initialized according to each element. We can see this each co-location point and you can see this we access the particular element here with this notation and finally assemble it to the Global Form. We can also look at this graphically, actually exactly at the point inside the stiffness matrix where the elements, the velocity changes and this can be seen here. You can see this by the changing color here, this is the central band of the stiffness matrix. This is the section here, the part where actually the low velocity zone starts which is why the values of the stiffness matrix changes. Now, let's do the simulation. The extrapolation scheme is exactly the same. You can see here in the graph that we highlighted in light blue the area in which we have a low velocity zone where we inject the source at the center. Now, as we hit the boundary, you can see this very nicely here, there is a reflection and transmission which can be expressed with reflection and transmission coefficients. They depend of course on the changes of the physical parameters. Now you can see actually similar to what we've seen in the finite difference simulation of a fault zone, that we get back and forth bouncing waves inside that physical domain. We also get actually reflections from the domain boundaries here at zero and 8,000 meters. Their free surface boundaries are perfectly reflecting and the simulation goes on. We can now see more and more phases appearing as transmissions and reflections from all the interfaces that we have already in the model. Now, a quite nice and elegant way of looking at this is actually this plot here, which is basically the whole wavefield shown in space-time. In this case, time actually starts from the top. So, the horizontal line is basically the physical domain in x from zero to eight kilometers, and from top to bottom time is increasing. Again, the amplitude scale is not important here. But you can see now nicely that from the central source point, the waves are propagating outwards and the slope of these wave fronts here in space-time of course depend on the velocities. So, you have this like the smaller the slope here to some extent, the faster the velocity. You can see these phases also they change polarity. That's also a famous effect. Here they change polarity and change polarity back, while they do not change polarity at the free surface only at these low velocity, high velocity boundaries. So, it's an elegant way of displaying the wave field here in space-time. So, this basically concludes our Python simulations. This was to some extent the most complicated simulation we've done, this was using the spectral element method on a heterogeneous velocity model. There is one interesting extension that we could use here that's a powerful aspect that we can apply for the spectral element method, that is making the element size vary. That's actually very easy to implement and I will leave this as an exercise and you can already imagine what that means, it basically means the Jacobian, which we have to use for the calculation of each elemental stiffness and mass matrix that will now become dependent of space. It's relatively straightforward to implement. But even for these low velocity zone implementation that also leads to an improvement in terms of computational efficiency.