diff --git a/example/2D-advection-diffusion.F90 b/example/2D-advection-diffusion.F90
index d197640..1255204 100644
--- a/example/2D-advection-diffusion.F90
+++ b/example/2D-advection-diffusion.F90
@@ -12,7 +12,7 @@ pure function scalar_field(x,y) result(gaussian)
     double precision, intent(in) :: x(:), y(:)
     double precision gaussian(size(x),size(y))
     double precision, parameter :: pi = acos(-1D0)
-    double precision, parameter :: x0 = 1D0, y0 = 1D0, sigma = pi/5
+    double precision, parameter :: x0 = -pi/2, y0 = -pi/2, sigma = pi/8
     do concurrent(integer :: j=1:size(y)) default(none) shared(x,y,gaussian)
       associate(r => sqrt((x-x0)**2 + (y(j)-y0)**2))
         gaussian(:,j) = exp(-(r**2)/(2*sigma**2))
@@ -20,21 +20,21 @@ pure function scalar_field(x,y) result(gaussian)
     end do
   end function
 
-  pure function stagnation_point_velocity(x,y) result(grad_phi)
+  pure function taylor_green_velocity(x,y) result(velocity)
     double precision, intent(in) :: x(:), y(:)
-    double precision grad_phi(size(x),size(y),space_dimension)
+    double precision velocity(size(x),size(y),space_dimension)
     do concurrent(integer :: i=1:size(x), j=1:size(y))
-      grad_phi(i,j,:) = [x(i), -y(j)]
+      velocity(i,j,:) = [10*sin(x(i))*cos(y(j)), -10*cos(x(i))*sin(y(j))]
     end do
   end function
 
 end module
 
 program advection_diffusion_2D
-  !! Solve the advection-diffusion equation for a passive scalar moving through a static
-  !! 2D velocity field.
+  !! Solve the advection-diffusion equation for a passive scalar moving through a
+  !! static velocity field define by 2D Taylor-Green vortices.
   use julienne_m, only : file_t
-  use fields_m, only : scalar_field, stagnation_point_velocity
+  use fields_m, only : scalar_field, taylor_green_velocity
   use formal_m, only : scalar_2D_t, vector_2D_t, scalar_2D_initializer_i, vector_2D_initializer_i
   implicit none
 
@@ -42,28 +42,37 @@ program advection_diffusion_2D
   procedure(scalar_2D_initializer_i), pointer :: scalar_2D_initializer
   procedure(vector_2D_initializer_i), pointer :: velocity_2D_initializer
   type(scalar_2D_t) s
+  double precision, parameter :: pi = acos(-1D0)
 
   scalar_2D_initializer => scalar_field
-  velocity_2D_initializer => stagnation_point_velocity
+  velocity_2D_initializer => taylor_green_velocity
 
-  s = scalar_2D_t(scalar_2D_initializer, order=4, cells=[20,20], x_min=[-2D0,-2D0], x_max=[2D0,2D0])
+  s = scalar_2D_t(scalar_2D_initializer, order=4, cells=[51,51], x_min=[-pi,-pi], x_max=[pi,pi])
 
   associate(v => vector_2D_t(velocity_2D_initializer, mold=s))
 
+    associate( &
+       scalar_file   => s%to_file("scalar") &
+      ,velocity_file => v%to_file("velocity") &
+    )
+      call   scalar_file%write_lines("example/scripts/scalar-initial.csv")
+      call velocity_file%write_lines("example/scripts/velocity.csv")
+    end associate
+
     advance_time: &
     block
-      double precision :: dt = 1D-6
-      associate(s_half => s + (dt/2) * d_dt(s, v))
-        s = s + dt * d_dt(s_half, v)
-      end associate
+      double precision :: dt = 1D-4
+      integer step
+
+      do step = 1, 500
+        associate(s_half => s + (dt/2) * d_dt(s, v))
+          s = s + dt * d_dt(s_half, v)
+        end associate
+      end do
     end block advance_time
 
-    associate( &
-       scalar_file   => s%to_file("scalar") &
-      ,velocity_file => v%to_file("vector") &
-    )
-      call   scalar_file%write_lines("example/scripts/scalar-adv-dif.csv")
-      call velocity_file%write_lines("example/scripts/velocity-adv-dif.csv")
+    associate(scalar_file => s%to_file("scalar"))
+      call scalar_file%write_lines("example/scripts/scalar-final.csv")
     end associate
 
   end associate
@@ -74,7 +83,7 @@ pure function d_dt(s, v) result(ds_dt)
     type(scalar_2D_t), intent(in) :: s
     type(vector_2D_t), intent(in) :: v
     type(scalar_2D_t) ds_dt
-    double precision, parameter :: D = 1D0
+    double precision, parameter :: D = 0.2D0
     ds_dt = .div. (D * .grad. s) - .div. (v * s)
   end function
 
diff --git a/example/scripts/2D-scalar-field.gnuplot b/example/scripts/2D-scalar-field.gnuplot
index 1c7b965..ed8cab0 100644
--- a/example/scripts/2D-scalar-field.gnuplot
+++ b/example/scripts/2D-scalar-field.gnuplot
@@ -12,13 +12,20 @@ datafile = base_name . ".csv"
 
 set datafile separator ","
 
-# --- 1. Read column headers from line 1 ---
-xlabel = "" ; ylabel = "" ; zlabel = ""
-set table $Dummy
-    plot datafile every ::0::0 \
-        using (xlabel=strcol(1), ylabel=strcol(2), zlabel=strcol(3), 0):(0) \
-        with table
-unset table
+# --- 1. Read column headers from line 1 directly via the shell ------------
+#     The data is split into blank-line-separated x-slices (needed so pm3d
+#     draws the surface correctly). Because of that, gnuplot's own
+#     "every ::0::0" doesn't just grab line 1: with no block restriction it
+#     samples the first point of *every* slice, and the assignments in the
+#     "using" clause just get overwritten slice by slice -- what's left at
+#     the end is whatever the last slice's first point happened to be,
+#     which is exactly the garbled numeric title ("0.180...E-34(-3.14...,
+#     3.14...)") you were seeing. Reading the header straight off disk with
+#     the shell sidesteps that entirely.
+get_field(n) = system("head -n 1 " . datafile . " | awk -F',' -v n=" . n . " '{v=$n; gsub(/^[ \\t]+|[ \\t]+$/,\"\",v); print v}'")
+xlabel = get_field(1)
+ylabel = get_field(2)
+zlabel = get_field(3)
 
 # --- 2. Plot ---
 set title  zlabel . "(" . xlabel . ", " . ylabel . ")"
@@ -32,6 +39,9 @@ set ticslevel 0 ; set key off
 set terminal gif size 800,600
 set output base_name . ".gif"
 
-splot datafile every ::1 using 1:2:3 with pm3d title ""
+# Header is skipped via a shell "tail" pipe rather than "every ::1", so the
+# blank lines separating x-slices are preserved (pm3d still needs them) and
+# no per-slice point gets silently dropped the way "every ::1" was doing.
+splot "< tail -n +2 " . datafile . "" using 1:2:3 with pm3d title ""
 
 set output    # flush and close the file