diff --git a/docs/source/api/optical_elements/axiparabola.rst b/docs/source/api/optical_elements/axiparabola.rst
new file mode 100644
index 00000000..a7e171a2
--- /dev/null
+++ b/docs/source/api/optical_elements/axiparabola.rst
@@ -0,0 +1,5 @@
+Axiparabola
+===========
+
+.. autoclass:: lasy.optical_elements.Axiparabola
+    :members:
diff --git a/docs/source/api/optical_elements/index.rst b/docs/source/api/optical_elements/index.rst
index ec30451f..770d187e 100644
--- a/docs/source/api/optical_elements/index.rst
+++ b/docs/source/api/optical_elements/index.rst
@@ -6,3 +6,4 @@ Optical elements
 
    parabolic_mirror
    polynomial_spectral_phase
+   axiparabola
diff --git a/docs/source/tutorials/axiparabola.ipynb b/docs/source/tutorials/axiparabola.ipynb
new file mode 100644
index 00000000..21c85dd9
--- /dev/null
+++ b/docs/source/tutorials/axiparabola.ipynb
@@ -0,0 +1,241 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "996ae31d-192b-4988-bd6b-649c833ae967",
+   "metadata": {},
+   "source": [
+    "# Initializing a flying-focus laser from an axiparabola"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f25efcc9-1ed0-4f9a-8c62-da36f56ba02f",
+   "metadata": {},
+   "source": [
+    "In this example, we generate a \"flying-focus\" laser from an axiparabola. This is done by sending a super-Gaussian laser (in the near-field) onto an axiparabola and propagating it to the far field."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7f076564-8330-4af0-ad92-708e1825596f",
+   "metadata": {},
+   "source": [
+    "## Generate a super-Gaussian laser"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "eb5e1364-1b92-42de-9518-dedf62be4544",
+   "metadata": {},
+   "source": [
+    "Define the physical profile, as combination of a longitudinal and transverse profile."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ed0a409c-3c02-4fc4-b41b-ef41d0fb5a3a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from lasy.laser import Laser\n",
+    "from lasy.profiles.gaussian_profile import CombinedLongitudinalTransverseProfile\n",
+    "from lasy.profiles.longitudinal import GaussianLongitudinalProfile\n",
+    "from lasy.profiles.transverse import SuperGaussianTransverseProfile\n",
+    "\n",
+    "wavelength     = 800e-9  # Laser wavelength in meters\n",
+    "polarization   = (1,0)   # Linearly polarized in the x direction\n",
+    "energy         = 1.5     # Energy of the laser pulse in joules\n",
+    "spot_size      = 1e-3    # Spot size in the near-field: millimeter-scale\n",
+    "pulse_duration = 30e-15  # Pulse duration of the laser in seconds\n",
+    "t_peak         = 0.0     # Location of the peak of the laser pulse in time\n",
+    "\n",
+    "laser_profile = CombinedLongitudinalTransverseProfile(\n",
+    "    wavelength, polarization, energy,\n",
+    "    GaussianLongitudinalProfile(wavelength, pulse_duration, t_peak),\n",
+    "    SuperGaussianTransverseProfile(spot_size, n_order=16)\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c303d3b6-2422-445b-a332-bc3f48388612",
+   "metadata": {},
+   "source": [
+    "Define the grid on which this profile is evaluated. \n",
+    "\n",
+    "**The grid needs to be wide enough to contain the millimeter-scale spot size, but also fine enough to resolve the micron-scale laser wavelength.**"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "84cb6fd2-2f44-4cc9-a5c1-6dd1e0d72c67",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dimensions     = 'rt'                              # Use cylindrical geometry\n",
+    "lo             = (0,-2.5*pulse_duration)           # Lower bounds of the simulation box\n",
+    "hi             = (1.1*spot_size,2.5*pulse_duration)  # Upper bounds of the simulation box\n",
+    "num_points     = (3000, 30)                         # Number of points in each dimension\n",
+    "\n",
+    "laser = Laser(dimensions,lo,hi,num_points,laser_profile)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bb8dc634-b3c4-4448-834b-fe9acb6c9645",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "laser.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1f588477-90c0-45c3-a47c-2b57733b462f",
+   "metadata": {},
+   "source": [
+    "## Propagate the laser through the axiparabola, and to the far field.\n",
+    "\n",
+    "First, define the parameters of the axiparabola."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e3eae91f-edb7-44ac-9b92-5c1a1328a4b3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from lasy.optical_elements import Axiparabola\n",
+    "f0 = 3e-2     # Focal distance\n",
+    "delta = 1.5e-2  # Focal range\n",
+    "R = spot_size # Radius\n",
+    "axiparabola = Axiparabola( f0, delta, R )"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "292c6e63-0480-4fb9-b1ad-6b679a3472d0",
+   "metadata": {},
+   "source": [
+    "Apply the effect of the axiparabola, and then propagate the laser for a distance `z=f0` (beginning of the focal range)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3158e51d-e331-438d-90c1-5f8c63bf9841",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "laser.apply_optics( axiparabola )\n",
+    "laser.propagate( f0 )"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fa581cef-fc3d-414b-9f19-df8545d3b62f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "laser.show()\n",
+    "plt.ylim(-0.25e-3, 0.25e-3)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a5582bfe-ce2d-4aec-a757-03a48d4f16f2",
+   "metadata": {},
+   "source": [
+    "At this point, the laser can be saved to file, and used e.g. as input to a PIC simulation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9b8c068a-715b-400f-b6ed-7711a60c0f1f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "laser.write_to_file('flying_focus', 'h5')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "192cabcd-b96b-4bff-80b0-2c59d5a780cb",
+   "metadata": {},
+   "source": [
+    "## Check that the electric field on axis remains high over many Rayleigh ranges\n",
+    "\n",
+    "An axiparabola can maintain a high laser field over a long distance (larger than the Rayleigh length).\n",
+    "Here, we can check that the laser field remains high over several Rayleigh length."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ef0d6f2e-3aac-47a0-baef-39b42d4c5749",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import math\n",
+    "ZR = math.pi*wavelength*f0**2/spot_size**2\n",
+    "print('Rayleigh length: %.f mm' %(1.e3*ZR))\n",
+    "print('Focal range: %.f mm' %(1.e3*delta))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e2113f14-3ec7-4392-8d91-96f6e28dcc71",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "laser.propagate(2*ZR)\n",
+    "laser.show()\n",
+    "plt.ylim(-0.25e-3, 0.25e-3)\n",
+    "plt.title('Laser field after 2 Rayleigh range')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3cb4bd3b-ad8d-471f-a2c0-1257c1b3bd39",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "laser.propagate(2*ZR)\n",
+    "laser.show()\n",
+    "plt.ylim(-0.25e-3, 0.25e-3)\n",
+    "plt.title('Laser field after 4 Rayleigh range')"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/docs/source/tutorials/index.rst b/docs/source/tutorials/index.rst
index ed7c4913..d03510c5 100644
--- a/docs/source/tutorials/index.rst
+++ b/docs/source/tutorials/index.rst
@@ -5,3 +5,4 @@ Tutorials
    :maxdepth: 1
 
    gaussian_laser.ipynb
+   axiparabola.ipynb
diff --git a/lasy/optical_elements/__init__.py b/lasy/optical_elements/__init__.py
index 3a622c84..69b5bf0d 100644
--- a/lasy/optical_elements/__init__.py
+++ b/lasy/optical_elements/__init__.py
@@ -1,4 +1,5 @@
+from .axiparabola import Axiparabola
 from .parabolic_mirror import ParabolicMirror
 from .polynomial_spectral_phase import PolynomialSpectralPhase
 
-__all__ = ["ParabolicMirror", "PolynomialSpectralPhase"]
+__all__ = ["ParabolicMirror", "PolynomialSpectralPhase", "Axiparabola"]
diff --git a/lasy/optical_elements/axiparabola.py b/lasy/optical_elements/axiparabola.py
new file mode 100644
index 00000000..adeb22de
--- /dev/null
+++ b/lasy/optical_elements/axiparabola.py
@@ -0,0 +1,71 @@
+import numpy as np
+from scipy.constants import c
+
+from .optical_element import OpticalElement
+
+
+class Axiparabola(OpticalElement):
+    r"""
+    Class that represents an axiparabola.
+
+    The rays that impinge on the axiparabola at different radii are focused
+    to different positions on the axis (resulting in an extended "focal range").
+
+    For more details, see S. Smartsev et al, "Axiparabola: a long-focal-depth,
+    high-resolution mirror for broadband high-intensity lasers", Optics Letters 44, 14 (2019)
+
+    Parameters
+    ----------
+    f0: float (in meter)
+        The focal distance, i.e. the distance, from the axiparabola,
+        where the focal range starts.
+
+    delta: float (in meter)
+        The length of the focal range.
+
+    R: float (in meter)
+        The radius of the axiparabola. Rays coming from ``r=0`` focus
+        at ``z=f0`` ; rays coming from ``r=R`` focus at ``z=f0+delta``
+    """
+
+    def __init__(self, f0, delta, R):
+        self.f0 = f0
+        self.delta = delta
+        self.R = R
+
+    def amplitude_multiplier(self, x, y, omega, omega0):
+        """
+        Return the amplitude multiplier.
+
+        Parameters
+        ----------
+        x, y, omega : ndarrays of floats
+            Define points on which to evaluate the multiplier.
+            These arrays need to all have the same shape.
+        omega0 : float (in rad/s)
+            Central angular frequency, as used for the definition
+            of the laser envelope.
+
+        Returns
+        -------
+        multiplier : ndarray of complex numbers
+            Contains the value of the multiplier at the specified points.
+            This array has the same shape as the array omega.
+        """
+        # Implement Eq. 4 in Smatsev et al.
+        r2 = x**2 + y**2
+        sag = (
+            (1.0 / (4 * self.f0)) * r2
+            - (self.delta / (8 * self.f0**2 * self.R**2)) * r2**2
+            + self.delta
+            * (self.R**2 + 8 * self.f0 * self.delta)
+            / (96 * self.f0**4 * self.R**4)
+            * r2**3
+        )
+
+        # Calculate phase shift
+        T = np.exp(-2j * (omega / c) * sag)
+        # Remove intensity beyond R
+        T[x**2 + y**2 > self.R**2] = 0
+
+        return T
diff --git a/tests/test_axiparabola.py b/tests/test_axiparabola.py
new file mode 100644
index 00000000..21678d97
--- /dev/null
+++ b/tests/test_axiparabola.py
@@ -0,0 +1,75 @@
+"""
+This test checks the implementation of the axiparabola
+by initializing a super-Gaussian pulse in the near field, and
+propagating it to the middle of the focal range. It then
+checks that the field amplitude remains high over the focal range.
+"""
+
+import numpy as np
+
+from lasy.laser import Laser
+from lasy.optical_elements import Axiparabola
+from lasy.profiles.gaussian_profile import CombinedLongitudinalTransverseProfile
+from lasy.profiles.longitudinal import GaussianLongitudinalProfile
+from lasy.profiles.transverse import SuperGaussianTransverseProfile
+
+
+def test_axiparabola():
+
+    # Define the laser profile
+    wavelength = 800e-9  # Laser wavelength in meters
+    polarization = (1, 0)  # Linearly polarized in the x direction
+    energy = 1.5  # Energy of the laser pulse in joules
+    spot_size = 1e-3  # Spot size in the near-field: millimeter-scale
+    pulse_duration = 30e-15  # Pulse duration of the laser in seconds
+    t_peak = 0.0  # Location of the peak of the laser pulse in time
+    laser_profile = CombinedLongitudinalTransverseProfile(
+        wavelength,
+        polarization,
+        energy,
+        GaussianLongitudinalProfile(wavelength, pulse_duration, t_peak),
+        SuperGaussianTransverseProfile(spot_size, n_order=16),
+    )
+
+    # Define the laser on a grid
+    dimensions = "rt"  # Use cylindrical geometry
+    lo = (0, -2.5 * pulse_duration)  # Lower bounds of the simulation box
+    hi = (1.1 * spot_size, 2.5 * pulse_duration)  # Upper bounds of the simulation box
+    num_points = (3000, 30)  # Number of points in each dimension
+    laser = Laser(dimensions, lo, hi, num_points, laser_profile)
+
+    # Define the parameters of the axiparabola
+    f0 = 3e-2  # Focal distance
+    delta = 1.5e-2  # Focal range
+    R = spot_size  # Radius
+    axiparabola = Axiparabola(f0, delta, R)
+
+    # Apply axiparabola and propagate to the middle of the focal range
+    laser.apply_optics(axiparabola)
+    laser.propagate(f0 + delta / 2)
+    E_middle = abs(laser.grid.get_temporal_field()).max()
+
+    # Compute the equivalent Rayleigh length
+    import math
+
+    ZR = math.pi * wavelength * f0**2 / spot_size**2
+    # Check that the focal range is much longer than the Rayleigh length
+    assert delta > 5 * ZR
+
+    laser.propagate(
+        -2 * ZR
+    )  # Propagate to two Rayleigh lengths before the middle of the focal range
+    E_before = abs(laser.grid.get_temporal_field()).max()
+    # For a Gaussian beam, the field amplitude should be reduced
+    # by a factor of sqrt(5) after two Rayleigh lengths
+    # Here, we check that the field remains significantly higher (1.5 times higher)
+    assert E_before > E_middle / np.sqrt(5) * 1.5
+
+    laser.propagate(
+        4 * ZR
+    )  # Propagate to two Rayleigh lengths after the middle of the focal range
+    E_after = abs(laser.grid.get_temporal_field()).max()
+    # For a Gaussian beam, the field amplitude should be reduced
+    # by a factor of sqrt(5) after two Rayleigh lengths
+    # Here, we check that the field remains significantly higher (1.5 times higher)
+    assert E_after > E_after / np.sqrt(5) * 1.5