diff --git a/docs/make.jl b/docs/make.jl
index dd28a33..6de865b 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -13,7 +13,12 @@ makedocs(;
     ),
     pages=[
         "Home" => "index.md",
-        "Examples" => "man/examples.md",
+        "Theory" => [
+            "Simulations" => "man/simulations.md",
+            "Darboux Transformations" =>  "man/dt.md"],
+        "Examples" => [
+            "Simulation Examples" => "man/sim_examples.md",
+            "DT Examples" =>  "man/dt_examples.md"],
         "Visualization" => "man/visualization.md",
         "Library" => [
             "Public API" => "public.md",
diff --git a/docs/src/index.md b/docs/src/index.md
index 559541f..c6e3ef6 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -15,20 +15,21 @@ The following features are currently available:
 - Computing the integrals of motion (energy, momentum, and particle number) and their errors.
 - Computing the Darboux Transformation to study complicated analytical solutions. We currently support the breather and soliton seeds for extended nonlinear Schrodinger equations of order up to 5 (including cubic NLS, Hirota, LPD, Quintic, and arbitrary combinations thereof). We also support the `cn` and `dn` seeds for the cubic NLS.
 - Easy [Visualization](@ref) through `Plots.jl` recipes.
-- Very simple API that allows one to compute very complicated solutions via only a few lines of code. Some examples can be found on the [Examples](@ref) page.
+- Very simple API that allows one to compute very complicated solutions via only a few lines of code.
 
 
+## Examples
 ```@contents
 Pages = [
-    "man/theory.md",
+    "man/dt_examples.md",
 ]
-Depth = 1
+Depth = 2
 ```
-
-## Library Outline
-
 ```@contents
-Pages = ["public.md"]
+Pages = [
+    "man/sim_examples.md",
+]
+Depth = 2
 ```
 
 ### [Index](@id main-index)
diff --git a/docs/src/man/dt.md b/docs/src/man/dt.md
new file mode 100644
index 0000000..a4db973
--- /dev/null
+++ b/docs/src/man/dt.md
@@ -0,0 +1,25 @@
+# Darboux Transformations
+
+Please see the paper to review the theory of Darboux Trasformations, and the [Darboux Transformation Examples](@ref) page for examples.
+
+## Supported Equations
+
+The most general equation supported by the Darboux transformation in `NonlinearSchrodinger.jl` is of the following form:
+
+```math
+i{\psi _x} + S[\psi (x,t)] - i\alpha H[\psi (x,t)] + \gamma P[\psi (x,t)] - i\delta Q[\psi (x,t)] = 0,
+```
+
+where
+```math
+\begin{aligned}
+S[\psi (x,t)] &= \frac{1}{2}{\psi _{tt}} + {\left| \psi  \right|^2}\psi, \\
+H[\psi (x,t)] &= {\psi _{ttt}} + 6{\left| \psi  \right|^2}{\psi _t}, \\
+P[\psi (x,t)] &= {\psi _{tttt}} + 8{\left| \psi  \right|^2}{\psi _{tt}} + 6{\left| \psi  \right|^4}\psi + 4{\left| {{\psi _t}} \right|^2}\psi + 6{\psi _t}^2{\psi ^*} + 2{\psi ^2}\psi _{tt}^*, \\
+Q[\psi (x,t)] &= {\psi _{ttttt}} + 10{\left| \psi  \right|^2}{\psi _{ttt}} + 30{\left| \psi  \right|^4}{\psi _t} + 10\psi {\psi _t}\psi _{tt}^* + 10\psi \psi _t^*{\psi _{tt}} + 20{\psi ^*}{\psi _t}{\psi _{tt}} + 10\psi _t^2\psi _t^*.
+\end{aligned}
+```
+
+Special cases include the cubic nonlinear Schrodinger equation (``\alpha = \gamma = \delta = 0``), the Hirota equation (``\alpha \neq 0, \gamma = \delta = 0``) the Lakshmanan-Porsezian-Daniel (LPD) equation (``\gamma \neq 0, \alpha = \delta = 0``) and the Quintic nonlinear Schrodinger equation (``\delta \neq 0, \alpha = \gamma = 0``).
+
+For this generalized NLSE, we support both the breather and soliton seeds. Additionally, for the cubic NLSE, we support the ``cn`` and ``dn`` seeds. Support for these seeds will be added at some point in the future for the generalized NLSE.
\ No newline at end of file
diff --git a/docs/src/man/dt_examples.md b/docs/src/man/dt_examples.md
new file mode 100644
index 0000000..716d8a1
--- /dev/null
+++ b/docs/src/man/dt_examples.md
@@ -0,0 +1,103 @@
+# Darboux Transformation Examples
+
+## Example 1: 7 Soliton Collision
+```@setup 1
+using NonlinearSchrodinger
+using Plots
+using LaTeXStrings
+```
+
+```@example 1
+xᵣ = -10=>10
+T = 20
+seed = "0"
+box = Box(xᵣ, T, Nₓ=500, Nₜ = 500)
+λ = [-0.45 + 0.775im, -0.35 + 0.8im, -0.25 + 0.825im, 0.85im, 0.25 + 0.875im, 0.35 + 0.9im, 0.45 + 0.925im]
+xₛ = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+tₛ = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+
+calc = Calc(λ, tₛ, xₛ, seed, box) 
+
+solve!(calc)
+heatmap(calc) 
+savefig("example1.svg") #hide
+```
+![](example1.svg)
+
+## Example 2: Fifth Order Maximal Intensity Breather
+```@setup 2
+using NonlinearSchrodinger
+using Plots
+using LaTeXStrings
+```
+
+```@example 2
+xᵣ = -10=>10
+λ₁ = 0.98im
+λ, T, Ω = params(λ = λ₁)
+box = Box(xᵣ, T, Nₓ=500, Nₜ = 500, n_periods = 3)
+
+λ = λ_maximal(λ₁, 5) # array of 5 eigenvalues
+xₛ = [0.0, 0.0, 0.0, 0.0, 0.0]
+tₛ = [0.0, 0.0, 0.0, 0.0, 0.0]
+
+seed = "exp"
+calc = Calc(λ, tₛ, xₛ, seed, box) 
+
+solve!(calc)
+surface(calc) 
+savefig("example2.svg") #hide
+```
+![](example2.svg)
+
+## Example 3: 3 Soliton Collision on a ``cn`` background
+```@setup 3
+using NonlinearSchrodinger
+using Plots
+using LaTeXStrings
+```
+
+```@example 3
+xᵣ = -10=>10
+T = 20
+
+box = Box(xᵣ, T, Nₓ=500, Nₜ = 500, n_periods = 1)
+λ = [-0.3+0.85im, 0.9im, 0.3+0.95im]
+xₛ = [0.0, 0.0, 0.0]
+tₛ = [0.0, 0.0, 0.0]
+
+seed = "cn"
+calc = Calc(λ, tₛ, xₛ, seed, box, m = 0.5) 
+
+solve!(calc)
+surface(calc) 
+savefig("example3.svg") #hide
+```
+![](example3.svg)
+
+## Example 4: First Order Breather matched to a ``dn`` Background
+```@setup 4 
+using NonlinearSchrodinger
+using Plots
+using LaTeXStrings
+```
+
+```@example 4
+xᵣ = -10=>10
+m = 2/5
+λ = λ_given_m(m, q=4)
+λ, T, Ω = params(λ = λ, m=m)
+box = Box(xᵣ, T, Nₓ=500, Nₜ = 500, n_periods = 3)
+
+λ = λ_maximal(λ, 1, m=m)
+xₛ = [0.0]
+tₛ = [0.0]
+
+seed = "dn"
+calc = Calc(λ, tₛ, xₛ, seed, box, m=m) 
+
+solve!(calc)
+heatmap(calc) 
+savefig("example4.svg") #hide
+```
+![](example4.svg)
\ No newline at end of file
diff --git a/docs/src/man/examples.md b/docs/src/man/examples.md
deleted file mode 100644
index 19f1093..0000000
--- a/docs/src/man/examples.md
+++ /dev/null
@@ -1,214 +0,0 @@
-# Examples
-
-## Example 1: Cosine Wave Initial Condition
-```@setup 1
-using NonlinearSchrodinger
-using Plots
-using LaTeXStrings
-```
-
-```@example 1
-λ, T, Ω = params(λ = 0.8im)
-
-xᵣ = 0=>20
-box = Box(xᵣ, T, dx=5e-3, Nₜ = 256, n_periods = 1)
-
-coeff = [1e-4]
-ψ₀, A₀ = ψ₀_periodic(coeff, box, Ω)
-
-sim = Sim(λ, box, ψ₀, T4A_TJ!)
-
-solve!(sim)
-compute_IoM!(sim)
-surface(sim)
-savefig("example1_psi.svg") #hide
-heatmap(sim, :ψ̃)
-savefig("example1_psi_tilde.svg") #hide
-plot(sim, :IoM)
-savefig("example1_IoM.svg") #hide
-```
-![](example1_psi.svg)
-![](example1_psi_tilde.svg)
-![](example1_IoM.svg)
-
-## Example 2: Soliton Initial Condition
-```@setup 2
-using NonlinearSchrodinger
-using Plots
-using LaTeXStrings
-```
-
-```@example 2
-λ = 0.75im
-
-T = 20
-xᵣ = 0=>20
-box = Box(xᵣ, T, dx=1e-2, Nₜ = 256, n_periods = 1)
-
-ψ₀ = Array{Complex{Float64}}(undef, box.Nₜ)
-ψ₀ .= 2*imag(λ)./cosh.(2*imag(λ).*box.t)
-
-sim = Sim(λ, box, ψ₀, T4A_TJ!)
-
-solve!(sim)
-surface(sim)
-savefig("example2.svg") # hide
-```
-![](example2.svg)
-
-## Example 3: Nonlinear Talbot Carpet (Pruning)
-```@setup 3
-using NonlinearSchrodinger
-using Plots
-using LaTeXStrings
-```
-
-```@example 3
-λ, T, Ω = params(a = 0.36)
-
-xᵣ = 0=>60
-box = Box(xᵣ, T, dx=1e-2, Nₜ = 512, n_periods = 5)
-
-coeff = [(2.7 + 4.6im)*1e-2]
-ψ₀, A₀ = ψ₀_periodic(coeff, box, Ω)
-
-sim = Sim(λ, box, ψ₀, T4A_TJ!, β=10.0)
-
-solve!(sim)
-heatmap(sim)
-savefig("example3_psi.svg") #hide
-plot(sim, :ψ̃)
-savefig("example3_psi_tilde.svg") #hide
-```
-![](example3_psi.svg)
-![](example3_psi_tilde.svg)
-
-## Example 4: 7 Soliton Collision
-```@setup 4
-using NonlinearSchrodinger
-using Plots
-using LaTeXStrings
-```
-
-```@example 4
-xᵣ = -10=>10
-T = 20
-seed = "0"
-box = Box(xᵣ, T, Nₓ=500, Nₜ = 500)
-λ = [-0.45 + 0.775im, -0.35 + 0.8im, -0.25 + 0.825im, 0.85im, 0.25 + 0.875im, 0.35 + 0.9im, 0.45 + 0.925im]
-xₛ = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-tₛ = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-
-calc = Calc(λ, tₛ, xₛ, seed, box) 
-
-solve!(calc)
-heatmap(calc) 
-savefig("example4.svg") #hide
-```
-![](example4.svg)
-
-## Example 5: Fifth Order Maximal Intensity Breather
-```@setup 5
-using NonlinearSchrodinger
-using Plots
-using LaTeXStrings
-```
-
-```@example 5
-xᵣ = -10=>10
-λ₁ = 0.98im
-λ, T, Ω = params(λ = λ₁)
-box = Box(xᵣ, T, Nₓ=500, Nₜ = 500, n_periods = 3)
-
-λ = λ_maximal(λ₁, 5) # array of 5 eigenvalues
-xₛ = [0.0, 0.0, 0.0, 0.0, 0.0]
-tₛ = [0.0, 0.0, 0.0, 0.0, 0.0]
-
-seed = "exp"
-calc = Calc(λ, tₛ, xₛ, seed, box) 
-
-solve!(calc)
-surface(calc) 
-savefig("example5.svg") #hide
-```
-![](example5.svg)
-
-## Example 6: 3 Soliton Collision on a ``cn`` background
-```@setup 6
-using NonlinearSchrodinger
-using Plots
-using LaTeXStrings
-```
-
-```@example 6
-xᵣ = -10=>10
-T = 20
-
-box = Box(xᵣ, T, Nₓ=500, Nₜ = 500, n_periods = 1)
-λ = [-0.3+0.85im, 0.9im, 0.3+0.95im]
-xₛ = [0.0, 0.0, 0.0]
-tₛ = [0.0, 0.0, 0.0]
-
-seed = "cn"
-calc = Calc(λ, tₛ, xₛ, seed, box, m = 0.5) 
-
-solve!(calc)
-surface(calc) 
-savefig("example6.svg") #hide
-```
-![](example6.svg)
-
-## Example 7: First Order Breather matched to a ``dn`` Background
-```@setup 7
-using NonlinearSchrodinger
-using Plots
-using LaTeXStrings
-```
-
-```@example 7
-xᵣ = -10=>10
-m = 2/5
-λ = λ_given_m(m, q=4)
-λ, T, Ω = params(λ = λ, m=m)
-box = Box(xᵣ, T, Nₓ=500, Nₜ = 500, n_periods = 3)
-
-λ = λ_maximal(λ, 1, m=m)
-xₛ = [0.0]
-tₛ = [0.0]
-
-seed = "dn"
-calc = Calc(λ, tₛ, xₛ, seed, box, m=m) 
-
-solve!(calc)
-heatmap(calc) 
-savefig("example7.svg") #hide
-```
-![](example7.svg)
-
-## Example 8: Darboux Transformation Initial Condition
-```@setup 8
-using NonlinearSchrodinger
-using Plots
-using LaTeXStrings
-```
-
-```@example 8
-λ₁ = 0.98im
-λ, T, Ω = params(λ = λ₁)
-
-xᵣ = 0=>20
-box = Box(xᵣ, T, dx=5e-3, Nₜ = 512, n_periods = 1)
-
-λ = λ_maximal(λ₁, 5) # array of 5 eigenvalues
-xₛ = [0.0, 0.0, 0.0, 0.0, 0.0]
-tₛ = [0.0, 0.0, 0.0, 0.0, 0.0]
-ψ₀ = ψ₀_DT(λ, tₛ, xₛ, -10, box)
-
-sim = Sim(λ₁, box, ψ₀, T4A_TJ!)
-
-solve!(sim)
-compute_IoM!(sim)
-heatmap(sim) 
-savefig("example8.svg") #hide
-```
-![](example8.svg)
\ No newline at end of file
diff --git a/docs/src/man/integrators.md b/docs/src/man/integrators.md
index b9dcc1d..d8caa91 100644
--- a/docs/src/man/integrators.md
+++ b/docs/src/man/integrators.md
@@ -1,5 +1,5 @@
 # Integrators
-This page shows the availabe integrators.
+This page shows the availabe integrators. Integrators prefixed with `_H` are used for the Hirota equation, and those with `_SS` are for the Sasa-Satsuma equation. Note that both these equations are not fully tested and should be used with utmost care due to the periodic boundary conditions.
 
 ```@index
 Pages = ["integrators.md"]
diff --git a/docs/src/man/sim_examples.md b/docs/src/man/sim_examples.md
new file mode 100644
index 0000000..f13a616
--- /dev/null
+++ b/docs/src/man/sim_examples.md
@@ -0,0 +1,112 @@
+# Simulation Examples
+
+## Example 1: Cosine Wave Initial Condition
+```@setup 1
+using NonlinearSchrodinger
+using Plots
+using LaTeXStrings
+```
+
+```@example 1
+λ, T, Ω = params(λ = 0.8im)
+
+xᵣ = 0=>20
+box = Box(xᵣ, T, dx=5e-3, Nₜ = 256, n_periods = 1)
+
+coeff = [1e-4]
+ψ₀, A₀ = ψ₀_periodic(coeff, box, Ω)
+
+sim = Sim(λ, box, ψ₀, T4A_TJ!)
+
+solve!(sim)
+compute_IoM!(sim)
+surface(sim)
+savefig("example1_psi.svg") #hide
+heatmap(sim, :ψ̃)
+savefig("example1_psi_tilde.svg") #hide
+plot(sim, :IoM)
+savefig("example1_IoM.svg") #hide
+```
+![](example1_psi.svg)
+![](example1_psi_tilde.svg)
+![](example1_IoM.svg)
+
+## Example 2: Soliton Initial Condition
+```@setup 2
+using NonlinearSchrodinger
+using Plots
+using LaTeXStrings
+```
+
+```@example 2
+λ = 0.75im
+
+T = 20
+xᵣ = 0=>20
+box = Box(xᵣ, T, dx=1e-2, Nₜ = 256, n_periods = 1)
+
+ψ₀ = Array{Complex{Float64}}(undef, box.Nₜ)
+ψ₀ .= 2*imag(λ)./cosh.(2*imag(λ).*box.t)
+
+sim = Sim(λ, box, ψ₀, T4A_TJ!)
+
+solve!(sim)
+surface(sim)
+savefig("example2.svg") # hide
+```
+![](example2.svg)
+
+## Example 3: Nonlinear Talbot Carpet (Pruning)
+```@setup 3
+using NonlinearSchrodinger
+using Plots
+using LaTeXStrings
+```
+
+```@example 3
+λ, T, Ω = params(a = 0.36)
+
+xᵣ = 0=>60
+box = Box(xᵣ, T, dx=1e-2, Nₜ = 512, n_periods = 5)
+
+coeff = [(2.7 + 4.6im)*1e-2]
+ψ₀, A₀ = ψ₀_periodic(coeff, box, Ω)
+
+sim = Sim(λ, box, ψ₀, T4A_TJ!, β=10.0)
+
+solve!(sim)
+heatmap(sim)
+savefig("example3_psi.svg") #hide
+plot(sim, :ψ̃)
+savefig("example3_psi_tilde.svg") #hide
+```
+![](example3_psi.svg)
+![](example3_psi_tilde.svg)
+
+## Example 4: Darboux Transformation Initial Condition
+```@setup 4
+using NonlinearSchrodinger
+using Plots
+using LaTeXStrings
+```
+
+```@example 4
+λ₁ = 0.98im
+λ, T, Ω = params(λ = λ₁)
+
+xᵣ = 0=>20
+box = Box(xᵣ, T, dx=5e-3, Nₜ = 512, n_periods = 1)
+
+λ = λ_maximal(λ₁, 5) # array of 5 eigenvalues
+xₛ = [0.0, 0.0, 0.0, 0.0, 0.0]
+tₛ = [0.0, 0.0, 0.0, 0.0, 0.0]
+ψ₀ = ψ₀_DT(λ, tₛ, xₛ, -10, box)
+
+sim = Sim(λ₁, box, ψ₀, T4A_TJ!)
+
+solve!(sim)
+compute_IoM!(sim)
+heatmap(sim) 
+savefig("example4.svg") #hide
+```
+![](example4.svg)
\ No newline at end of file
diff --git a/docs/src/man/simulations.md b/docs/src/man/simulations.md
new file mode 100644
index 0000000..4fe4395
--- /dev/null
+++ b/docs/src/man/simulations.md
@@ -0,0 +1,3 @@
+# Simulations
+
+To review the theory of how simulations are done using the split-step fourier method, please review the paper (to be posted on the arXiv soon). For examples, see the [Simulation Examples](@ref) page. For the available algorithms, see the [Integrators](@ref) page.
\ No newline at end of file
diff --git a/docs/src/man/visualization.md b/docs/src/man/visualization.md
index cf78458..453ae6c 100644
--- a/docs/src/man/visualization.md
+++ b/docs/src/man/visualization.md
@@ -1,7 +1,7 @@
 # Visualization
 
 ## Recipes
-`NonlinearSchrodinger.jl` provides recipes for `Plots.jl` to make plotting easy. You can simply call any plotting function on a `::Sim` or `::Calc` object, as shown in the [Examples](@ref) page.
+`NonlinearSchrodinger.jl` provides recipes for `Plots.jl` to make plotting easy. You can simply call any plotting function on a `::Sim` or `::Calc` object, as shown in the examples pages.
 
 ## Visualizing Recursion
 
diff --git a/src/Utilities.jl b/src/Utilities.jl
index 635d3eb..262dd11 100644
--- a/src/Utilities.jl
+++ b/src/Utilities.jl
@@ -4,7 +4,7 @@
 Computes the integrals of motion of `obj.ψ` and saves them in respective fields of `obj`.
 `obj` can be a `::Sim` or `::Calc` object. The fields are:
 
-````
+```
 obj.KE  # Array containing the kinetic energy K(x)
 obj.PE  # Array containing the potential energy V(x)
 obj.E   # Array containing the energy H(x)
@@ -13,7 +13,7 @@ obj.P   # Array containing the momentum P(x)
 obj.dE  # Array containing the energy error δE(x)
 obj.dN  # Array containing the norm error δN(x)
 obj.dP  # Array containing the momentum error δP(x)
-````
+```
 
 Result can be plotted using `plot(obj, :IoM)`
 """
@@ -47,7 +47,7 @@ end
 """
     function params(; kwargs...)
 
-Computes parameters
+Computes parameters `λ, T, Ω` given either of `λ, T, Ω, a` and possibly `m` for dnoidal seed. Used only for `exp` and `dn` seeds.
 """
 function params(;m=0.0, kwargs...)
     if length(kwargs) != 1
@@ -150,6 +150,13 @@ function ψ₀_periodic(coeff::Array, box::Box, Ω; phase=0)
     return ψ₀, A0
 end #psi0_periodic
 
+"""
+    ψ₀_DT(λ, tₛ, xₛ, X₀, box; seed="exp", f = Dict(:α=> 0.0, :γ => 0.0, :δ=>0.0))
+
+Computes an initial condition using the DT characterized by `λ, xₛ, tₛ` and `seed` at x = `X₀` for the extended equation characterized by a dictionary of parameters `f` (currently only `α`) is supported. `box::Box` is the simulation box.
+
+See also: [`Calc`](@ref), [`Sim`](@ref)
+"""
 function ψ₀_DT(λ, tₛ, xₛ, X₀, box; seed="exp", f = Dict(:α=> 0.0, :γ => 0.0, :δ=>0.0))
     xᵣ = X₀=>X₀+1e-5
     T = abs(box.t[1]*2)
@@ -162,6 +169,13 @@ function ψ₀_DT(λ, tₛ, xₛ, X₀, box; seed="exp", f = Dict(:α=> 0.0, :γ
     return calc.ψ[:, 1]
 end
 
+"""
+    λ_maximal(λ₁, N; m = 0)
+
+Computed a maximal intensity set of `λ` of order `N` given `λ₁` and possibly `m` for dnoidal bakcground.
+
+See also: [`λ_given_m`](@ref)
+"""
 function λ_maximal(λ₁, N; m = 0)
     ν₁ = imag(λ₁)
     #ν_min = sqrt(1 - 1/N^2)
@@ -178,11 +192,23 @@ function λ_maximal(λ₁, N; m = 0)
     λ = sqrt.(G_n .+ sqrt.(G_n.^2 .- 64*m^2*ν₁^4))./(4*sqrt(2)*ν₁)*im
 end
 
+"""
+    λ_given_m(m; q = 2)
+
+Computed a `λ` that is matched to the dnoidal background given by `m` with an integer `q`. See paper for more details.
+
+See also: [`λ_maximal`](@ref)
+"""
 function λ_given_m(m; q = 2)
     F = π/(2*q*Elliptic.K(m))
     λ = 0.5 * sqrt(2 - 2*F^2 - m + 2*sqrt((F^2-1)*(F^2-1+m)))*im
 end
 
+"""
+    λ_given_f(f, ν)
+
+Computes `λ = v + i ν` such that `v` results in breather to soliton conversion for the extended NLS characterized by the dictionary of parameters `f`.
+"""
 function λ_given_f(f, ν)
     α = f[:α]
     γ = f[:γ]
@@ -204,6 +230,11 @@ function λ_given_f(f, ν)
     end
 end
 
+"""
+    PHF(calc::Calc)
+
+Computed the peak height formula for any `calc::Calc` objecti irrespective of seed.
+"""
 function PHF(calc::Calc)
     s = 2*sum(imag.(calc.λ))
     if calc.seed == "exp" || calc.seed == "dn"
diff --git a/src/kawahara.jl b/src/kawahara.jl
new file mode 100644
index 0000000..2830705
--- /dev/null
+++ b/src/kawahara.jl
@@ -0,0 +1,32 @@
+using Plots
+using FFTW
+using NonlinearSchrodinger
+
+function T1A_K(ψᵢ, b)
+    γ = 1
+    p = 1
+    q = 1
+
+    n = 1
+    μ = 0 
+
+    #ψ = ifft(im*b.ω.*fftshift(fft(ψᵢ)))
+    #ψₒ = exp.(-γ*ψᵢ.^(n-1).*b.dx.*ψ).*ψᵢ
+    #ψₒ = fftshift(fft(ψₒ))
+    #ψₒ = ifft(exp.(-b.dx*(im*b.ω*μ .+ im*b.ω.^3*p .+ im*b.ω.^5*q)).*ψₒ)
+
+    ψₒ = fftshift(fft(ψᵢ))
+    ψₒ = ifft(exp.(-b.dx*(im*b.ω*μ .+ im*b.ω.^3*p .+ im*b.ω.^5*q)).*ψₒ)
+    ψₒ = exp.(-γ*ψᵢ.^(n-1).*b.dx.*ifft(+im*b.ω.*fftshift(fft(ψₒ)))).*ψₒ
+end
+
+xᵣ = 0=>50
+box = Box(xᵣ, 200, dx=5e-2, Nₜ = 512, n_periods =3)
+
+ψ = zeros(box.Nₓ, box.Nₜ) .+ im
+
+ψ[1, :] = 105/169*sech.(1/2/sqrt(13)*box.t).^4
+for i = 1:box.Nₓ-1
+   ψ[i+1, :] .=  T1A_K(ψ[i,:], box)
+end
+heatmap(abs.(ψ))
\ No newline at end of file