diff --git a/Project.toml b/Project.toml
index 40f3650..7a8e6df 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,17 +6,20 @@ version = "0.1.0"
 [deps]
 ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
-MathTeXEngine = "0a4f8689-d25c-4efe-a92b-7142dfc1aa53"
-FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MathTeXEngine = "0a4f8689-d25c-4efe-a92b-7142dfc1aa53"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 SymbolicNumericIntegration = "78aadeae-fbc0-11eb-17b6-c7ec0477ba9e"
+SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
diff --git a/misc/FDCoeff.ipynb b/misc/FDCoeff.ipynb
new file mode 100644
index 0000000..fc45184
--- /dev/null
+++ b/misc/FDCoeff.ipynb
@@ -0,0 +1,129 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 49,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$$ \\begin{equation}\n",
+       "\\left[\n",
+       "\\begin{array}{c}\n",
+       "f \\\\\n",
+       "\\end{array}\n",
+       "\\right]\n",
+       "\\end{equation}\n",
+       " $$"
+      ],
+      "text/plain": [
+       "1-element Vector{Num}:\n",
+       " f"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "using Symbolics\n",
+    "\n",
+    "\n",
+    "@variables x₁ x₂ x₃\n",
+    "@variables ∂₁ ∂₂ ∂₃\n",
+    "@variables ΔX₁ ΔX₂ ΔX₃\n",
+    "#@variables f(..)\n",
+    "@variables f\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 56,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "q += p = f + ΔX₁*∂₁ + ΔX₂*∂₂ + ΔX₃*∂₂\n"
+     ]
+    },
+    {
+     "ename": "ErrorException",
+     "evalue": "syntax: invalid iteration specification",
+     "output_type": "error",
+     "traceback": [
+      "syntax: invalid iteration specification\n",
+      "\n",
+      "Stacktrace:\n",
+      " [1] top-level scope\n",
+      "   @ ~/Github/SeismicQ/misc/FDCoeff.ipynb:9"
+     ]
+    }
+   ],
+   "source": [
+    "q =f\n",
+    "ordre = 2 # bigger than 2\n",
+    "p = ∂₁ * ΔX₁ + ∂₂ * ΔX₂ + ∂₂ * ΔX₃\n",
+    "q += p\n",
+    "\n",
+    "\n",
+    "\n",
+    "for m=2,ordre\n",
+    "    @show p *= 1//m * (∂₁ * ΔX₁ + ∂₂ * ΔX₂ + ∂₂ * ΔX₃)\n",
+    "    @show q += p\n",
+    "\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 35,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$$ \\begin{equation}\n",
+       "2 f\n",
+       "\\end{equation}\n",
+       " $$"
+      ],
+      "text/plain": [
+       "2f"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "q = f+f"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Julia 1.9.2",
+   "language": "julia",
+   "name": "julia-1.9"
+  },
+  "language_info": {
+   "file_extension": ".jl",
+   "mimetype": "application/julia",
+   "name": "julia",
+   "version": "1.9.2"
+  },
+  "orig_nbformat": 4
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/misc/Test_Symbolics_for_OptimallyAccurateOperators.jl b/misc/Test_Symbolics_for_OptimallyAccurateOperators.jl
index 535aa2c..24382c6 100644
--- a/misc/Test_Symbolics_for_OptimallyAccurateOperators.jl
+++ b/misc/Test_Symbolics_for_OptimallyAccurateOperators.jl
@@ -7,4 +7,55 @@
 #
 # Titi told me that I should learn how to code with Symbolics so let's go (it's 16h on Friday though)
 #  
-#  Nobu
\ No newline at end of file
+#  Nobu
+
+#using Symbolics, SymbolicNumericIntegration
+#using DifferentialEquations, ModelingToolkit
+
+using Symbolics
+
+#@parameters x₁ x₂ x₃
+@variables ∂₁ ∂₂ ∂₃
+@variables ΔX₁ ΔX₂ ΔX
+@variables f
+@variables m
+
+# m-th order
+
+@show 1//factorial(m)*∂₁^
+for n=
+
+@parameters t σ ρ β
+@variables x(t) y(t) z(t)
+D = Differential(t)
+
+eqs = [D(D(x)) ~ σ * (y - x),
+    D(y) ~ x * (ρ - z) - y,
+    D(z) ~ x * y - β * z]
+
+@named sys = ODESystem(eqs)
+sys = structural_simplify(sys)
+
+u0 = [D(x) => 2.0,
+    x => 1.0,
+    y => 0.0,
+    z => 0.0]
+
+p = [σ => 28.0,
+    ρ => 10.0,
+    β => 8 / 3]
+
+tspan = (0.0, 100.0)
+prob = ODEProblem(sys, u0, tspan, p, jac = true)
+sol = solve(prob)
+using Plots
+plot(sol, idxs = (x, y))
+
+
+# Symbols
+
+
+
+
+
+