--- a/lattice2d.jl	Wed Jan 05 18:21:11 2022 +0100
+++ b/lattice2d.jl	Wed Jan 05 21:18:21 2022 +0100
@@ -2,6 +2,22 @@
 using DifferentialEquations
 import LinearAlgebra
 
+using SparseArrays
+
+theme(:juno)
+
+# Masses schema:
+#
+#  7 -- 8 -- 9
+#  |  X |  X |
+#  4 -- 5 -- 6
+#  |  X |  X |
+#  1 -- 2 -- 3
+# 
+# where 'X' is a cross connection.
+# We use Born-van Karman boundary conditions, i.e.
+# a "Pacman world" where the area behaves like it is on a torus.
+
 function neighbors(n::Int, i::Int)::Tuple
     # Returns neighbors in bottom/top/right/left
     a = (i + n) % n^2
@@ -31,13 +47,18 @@
 #
 # This assumes constant spring constant k for all springs, obviously.
 # In the vector notation, oscillator i has its x component at index 2i-1 and its y component at 2i.
+# Additional terms for diagonal springs are introduced. The spring constants for this are lower.
 
-function coupling_matrix(n::Int, k::Float64)::Matrix
+function coupling_matrix(n::Int, k::Float64; damping::Float64=0., full=false)
     N = n^2
-    A = zeros(4N, 4N)
+    if full
+        A = zeros(4N, 4N)
+    else
+        A = spzeros(4N, 4N)
+    end
 
     # Link velocity to displacement.
-    A[1:2N, 2N+1:4N] = LinearAlgebra.I(2N)
+    A[1:2N, 2N+1:end] = LinearAlgebra.I(2N)
 
     for i in 1:N
         (b,t,r,l) = neighbors(n, i)
@@ -49,6 +70,7 @@
         y = 2i
 
         # TO DO: Introduce non-linear coupling
+        # TO DO: Introduce damping
 
         A[2N+x, 2r-1] += k
         A[2N+x, 2l-1] += k
@@ -68,6 +90,9 @@
         end
         A[2N+x, x] -= 3/2 * k
         A[2N+y, y] -= 3/2 * k
+
+        # damping
+        A[2N+1:end, 2N+1:end] -= LinearAlgebra.I(2N) * damping
     end
     A
 end
@@ -96,34 +121,43 @@
 
     xs += dxy[1,:]
     ys += dxy[2,:]
-    scatter(xs .* a, ys .* a, xlim=(-a, (n)*a), ylim=(-2a, (n)*a), title="Positions $(i)")
+    scatter(xs .* a, ys .* a, xlim=(-a, n*a), ylim=(-a, n*a), title="Positions $(i)")
 end
 
 function f2d(du, u, A, t)
     LinearAlgebra.mul!(du, A, u)
 end
 
-n = 5
+# There will be n x n oscillators.
+n = 10
 
 # These are two different initial configurations:
-
-# All elements on two orthogonal sides are displaced initially.
+#  All elements on two orthogonal sides are displaced initially.
 initial_xy_displacements = vcat(
     [(1, i, 0., .2) for i in 1:n],
     [(i, 2, .2, .0) for i in 1:n])
-# One element is displaced diagonally.
-one_corner_displacement = [(1,1,0.2,0.2)]
+#  One element is displaced diagonally.
+initial_corner_displacements = [(1,1,0.3,0.3), (n, 1, 0.3, -0.3)]
 
-u0 = initial_state(n, initial_xy_displacements)
+u0 = initial_state(n, initial_corner_displacements)
+
+# Essential part: the coupling matrix!
+coupling = coupling_matrix(n, 1.; damping=0.0001)
 
-coupling = coupling_matrix(n, 1.)
+# ODE formalities
+problem = ODEProblem(f2d, u0, (0, 40), coupling)
+@time solution = solve(problem)
 
-# Look at a few eigenvibrations
-ED = LinearAlgebra.eigen(coupling[2n^2+1:end, 1:2n^2])
+@gif for (j, u) in enumerate(solution.u)
+    plot_state(n, u, j)
+end
 
+plot(solution, vars=[(0, 1), (0, 2)])
+
+# Plot animations for all eigenmodes. Careful: there are n^2 eigenmodes!
 function plot_all_eigenmodes(n=5)
     # Enumerate all eigenmodes
-    coupling = coupling_matrix(n, 1.)
+    coupling = coupling_matrix(n, 1.; full=true)
     ED = LinearAlgebra.eigen(coupling[2n^2+1:end, 1:2n^2])
     @show size(coupling)
 

--- a/linear_chain.jl	Wed Jan 05 18:21:11 2022 +0100
+++ b/linear_chain.jl	Wed Jan 05 21:18:21 2022 +0100
@@ -4,6 +4,7 @@
 using Plots
 theme(:juno)
 
+
 function n_coupled_masses(n::Int, ks::Vector, masses::Vector)
     # Prepares coupling matrix.
     if length(ks) != (n+1)