There are a few packages you can use for solving the PDE. Makie is great for creating animations. Hereâs an example, using FiniteVolumeMethod.jl for solving the PDE.

```
using CairoMakie, DelaunayTriangulation, FiniteVolumeMethod, OrdinaryDiffEq, LinearSolve
# Setup the geometry and initial condition
N = 500
L = 2.5
tri = triangulate_rectangle(-L, L, -L, L, N, N; single_boundary=true)
c, d, e, f, k, shift = 2.0, -1.0, 1.0, 0.5, 1.2, 10.0
Z = [(c^2 - (x / e - d)^2 - (y / f)^2)^2 + k * (c + d - x / e)^3 - shift for (x, y) in each_point(tri)]
@. Z = 1 - max(sign(Z), 0.0)
# Define the PDE
geo = FVMGeometry(tri)
bc = BoundaryConditions(geo, Returns(0.0), :Dirichlet)
final_time = 0.25
prob = FVMProblem(geo, bc; initial_condition=Z, diffusion_function=Returns(1.0), final_time=final_time)
# Solve
secs = 15.0
Ît = final_time / (24secs)
sol = solve(prob, TRBDF2(linsolve=KLUFactorization()), saveat=Ît)
# Plot
j = Observable(1)
fig = Figure(fontsize=33)
ax = Axis3(fig[1, 1],
xlabel=L"x", ylabel=L"y", zlabel=L"z",
title=lift(_j -> L"t = %$(rpad(round(sol.t[_j], digits = 5), 7, '0'))", j),
titlealign=:left)
surface!(ax, LinRange(-L, L, N), LinRange(-L, L, N), @lift(reshape(sol.u[$j], N, N)))
# If you just want a 2D plot, use Axis instead and tricontourf!(ax, tri, @lift(sol.u[$j]))
record(fig, "heat_animation.mp4", eachindex(sol); framerate=24) do _j
j[] = _j
end
```