using Plots
h = 0.05;
N = 200;
const g = 9;
const ℓ = 1.46;
θ = zeros(N+1)
v = zeros(N+1)
θ[1] = 1.0; # this is θ₀
v[1] = 0.5; # this is v₀
for k in 1:N
θ[k+1] = θ[k] + h * v[k]
v[k+1] = v[k] - h * ( g / ℓ ) * sin(θ[k+1])
end
plot(θ,v,lc=:red, lw=1)

The flow of this method preserve the area and I would like to investigate this property by using some pictures inside the plot like in the following

Take as starting positions a number of points on a circle and compute the iterates, F^{kn}, of these points, for a fixed n, and k=1, 2, …kmax. Plot the interpolating curve of the corresponding points after each k.

Thanks @empet ! The code work perfectly in this case, but when using explicit Euler like

using Plots
h = 0.05;
N = 200;
const g = 9;
const ℓ = 1.46;
θ = zeros(N+1)
v = zeros(N+1)
θ[1] = 1.0; # this is θ₀
v[1] = 0.5; # this is v₀
for k in 1:N
θ[k+1] = θ[k] + h * v[k]
v[k+1] = v[k] - h * ( g / ℓ ) * sin(θ[k])
end
plot(θ,v)

using Plots
function F()
h = 0.05;
N = 200;
g = 9;
ℓ = 1.46;
θ = zeros(N+1)
v = zeros(N+1)
θ[1] = 1.0; # this is θ₀
v[1] = 0.5; # this is v₀
for k in 1:N
θ[k+1] = θ[k] + h * v[k]
v[k+1] = v[k] - h * ( g / ℓ ) * sin(θ[k])
end
return hcat(θ, v)
end
composition(F, n) = ∘(ntuple(_ -> F, n)...)
function circle(t,k)
[1+(inv(1+h^2))^k*cos(t), (inv(1+h^2))^k*sin(t)]
end
t=range(0, 2pi, 35)
pc=circle.(t;)
pl= plot(getindex.(pc, 1), getindex.(pc, 2), lc=:red, lw=1,
legend=false, size=(450, 450), aspect_ratio = :equal)
for _ in 1:5:100
pc=composition(F, 10).(pc)
plot!(pl, getindex.(pc, 1), getindex.(pc, 2), lc=:red, lw=1)
plot!(x[:,1],x[:,2])
end
display(pl)

Your discrete system defined today is not the same as that you posted initially, when you formulated your question.
The last one, defined today,

θ[k+1] = θ[k] + h * v[k]
v[k+1] = v[k] - h * ( g / ℓ ) * sin(θ[k])

is not defined by an area preserving map, because the determinant of its Jacobian matrix is not 1 (one). That’s why your plot reveals an expansive effect, as the number of iterations increases.

Your initial system I have used to illustrate the area preserving of a disk under the action of the map F. was:

θ[k+1] = θ[k] + h * v[k]
v[k+1] = v[k] - h * ( g / ℓ ) * sin(θ[k+1])