diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
new file mode 100644
index 0000000..01bfab9
--- /dev/null
+++ b/.JuliaFormatter.toml
@@ -0,0 +1 @@
+# See https://domluna.github.io/JuliaFormatter.jl/stable/ for a list of options
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
new file mode 100644
index 0000000..e274804
--- /dev/null
+++ b/.github/workflows/CI.yml
@@ -0,0 +1,35 @@
+name: CI
+on:
+  push:
+    branches:
+      - master
+    tags: ['*']
+  pull_request:
+concurrency:
+  # Skip intermediate builds: always.
+  # Cancel intermediate builds: only if it is a pull request build.
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }}
+jobs:
+  test:
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - '1.9'
+          - 'nightly'
+        os:
+          - ubuntu-latest
+        arch:
+          - x64
+    steps:
+      - uses: actions/checkout@v3
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      - uses: julia-actions/cache@v1
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1
diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml
new file mode 100644
index 0000000..cba9134
--- /dev/null
+++ b/.github/workflows/CompatHelper.yml
@@ -0,0 +1,16 @@
+name: CompatHelper
+on:
+  schedule:
+    - cron: 0 0 * * *
+  workflow_dispatch:
+jobs:
+  CompatHelper:
+    runs-on: ubuntu-latest
+    steps:
+      - name: Pkg.add("CompatHelper")
+        run: julia -e 'using Pkg; Pkg.add("CompatHelper")'
+      - name: CompatHelper.main()
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          COMPATHELPER_PRIV: ${{ secrets.DOCUMENTER_KEY }}
+        run: julia -e 'using CompatHelper; CompatHelper.main()'
diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml
new file mode 100644
index 0000000..2bacdb8
--- /dev/null
+++ b/.github/workflows/TagBot.yml
@@ -0,0 +1,31 @@
+name: TagBot
+on:
+  issue_comment:
+    types:
+      - created
+  workflow_dispatch:
+    inputs:
+      lookback:
+        default: 3
+permissions:
+  actions: read
+  checks: read
+  contents: write
+  deployments: read
+  issues: read
+  discussions: read
+  packages: read
+  pages: read
+  pull-requests: read
+  repository-projects: read
+  security-events: read
+  statuses: read
+jobs:
+  TagBot:
+    if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot'
+    runs-on: ubuntu-latest
+    steps:
+      - uses: JuliaRegistries/TagBot@v1
+        with:
+          token: ${{ secrets.GITHUB_TOKEN }}
+          ssh: ${{ secrets.DOCUMENTER_KEY }}
diff --git a/.github/workflows/register.yml b/.github/workflows/register.yml
new file mode 100644
index 0000000..5b7cd3b
--- /dev/null
+++ b/.github/workflows/register.yml
@@ -0,0 +1,16 @@
+name: Register Package
+on:
+  workflow_dispatch:
+    inputs:
+      version:
+        description: Version to register or component to bump
+        required: true
+jobs:
+  register:
+    runs-on: ubuntu-latest
+    permissions:
+        contents: write
+    steps:
+      - uses: julia-actions/RegisterAction@latest
+        with:
+          token: ${{ secrets.GITHUB_TOKEN }}
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..5a16984
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,3 @@
+/Manifest.toml
+/docs/Manifest.toml
+/docs/build/
diff --git a/.travis.yml b/.travis.yml
new file mode 100644
index 0000000..c7561b6
--- /dev/null
+++ b/.travis.yml
@@ -0,0 +1,38 @@
+# Documentation: http://docs.travis-ci.com/user/languages/julia
+language: julia
+notifications:
+  email: false
+julia:
+  - 1.9
+  - nightly
+os:
+  - linux
+arch:
+  - x64
+cache:
+  directories:
+    - ~/.julia/artifacts
+jobs:
+  fast_finish: true
+  allow_failures:
+    - julia: nightly
+  include:
+    - stage: Documentation
+      julia: 1
+      script:
+      - |
+        julia --project=docs -e '
+          using Pkg
+          Pkg.develop(PackageSpec(path=pwd()))
+          Pkg.instantiate()'
+      - |
+        julia --project=docs -e '
+          using Pkg
+          Pkg.develop(PackageSpec(path=pwd()))
+          Pkg.instantiate()
+          using Documenter: DocMeta, doctest
+          using ImaginaryOptics
+          DocMeta.setdocmeta!(ImaginaryOptics, :DocTestSetup, :(using ImaginaryOptics); recursive=true)
+          doctest(ImaginaryOptics)
+          include("docs/make.jl")'
+      after_success: skip
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..639969e
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2023 Marco Menarini <menarini.marco@gmail.com>
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/Project.toml b/Project.toml
new file mode 100644
index 0000000..4588557
--- /dev/null
+++ b/Project.toml
@@ -0,0 +1,16 @@
+name = "ImaginaryOptics"
+uuid = "ed0a9f63-0a1c-4465-a5fd-c9e4b6c1d448"
+authors = ["Marco Menarini <menarini.marco@gmail.com>"]
+version = "1.0.0-DEV"
+
+[deps]
+Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
+
+[compat]
+julia = "1.9"
+
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test"]
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..2c9c6d3
--- /dev/null
+++ b/README.md
@@ -0,0 +1,7 @@
+# ImaginaryOptics
+
+[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://marcom.github.io/ImaginaryOptics.jl/stable/)
+[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://marcom.github.io/ImaginaryOptics.jl/dev/)
+[![Build Status](https://github.com/marcom/ImaginaryOptics.jl/actions/workflows/CI.yml/badge.svg?branch=master)](https://github.com/marcom/ImaginaryOptics.jl/actions/workflows/CI.yml?query=branch%3Amaster)
+[![Build Status](https://travis-ci.com/marcom/ImaginaryOptics.jl.svg?branch=master)](https://travis-ci.com/marcom/ImaginaryOptics.jl)
+[![PkgEval](https://JuliaCI.github.io/NanosoldierReports/pkgeval_badges/I/ImaginaryOptics.svg)](https://JuliaCI.github.io/NanosoldierReports/pkgeval_badges/I/ImaginaryOptics.html)
diff --git a/docs/Project.toml b/docs/Project.toml
new file mode 100644
index 0000000..5c372f1
--- /dev/null
+++ b/docs/Project.toml
@@ -0,0 +1,3 @@
+[deps]
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+ImaginaryOptics = "ed0a9f63-0a1c-4465-a5fd-c9e4b6c1d448"
diff --git a/docs/make.jl b/docs/make.jl
new file mode 100644
index 0000000..8b4b7c8
--- /dev/null
+++ b/docs/make.jl
@@ -0,0 +1,25 @@
+using ImaginaryOptics
+using Documenter
+
+DocMeta.setdocmeta!(ImaginaryOptics, :DocTestSetup, :(using ImaginaryOptics); recursive=true)
+
+makedocs(;
+    modules=[ImaginaryOptics],
+    authors="Marco Menarini <menarini.marco@gmail.com>",
+    repo="https://github.com/marcom/ImaginaryOptics.jl/blob/{commit}{path}#{line}",
+    sitename="ImaginaryOptics.jl",
+    format=Documenter.HTML(;
+        prettyurls=get(ENV, "CI", "false") == "true",
+        canonical="https://marcom.github.io/ImaginaryOptics.jl",
+        edit_link="master",
+        assets=String[],
+    ),
+    pages=[
+        "Home" => "index.md",
+    ],
+)
+
+deploydocs(;
+    repo="github.com/marcom/ImaginaryOptics.jl",
+    devbranch="master",
+)
diff --git a/docs/src/index.md b/docs/src/index.md
new file mode 100644
index 0000000..5dadd05
--- /dev/null
+++ b/docs/src/index.md
@@ -0,0 +1,14 @@
+```@meta
+CurrentModule = ImaginaryOptics
+```
+
+# ImaginaryOptics
+
+Documentation for [ImaginaryOptics](https://github.com/marcom/ImaginaryOptics.jl).
+
+```@index
+```
+
+```@autodocs
+Modules = [ImaginaryOptics]
+```
diff --git a/src/ImaginaryOptics.jl b/src/ImaginaryOptics.jl
new file mode 100644
index 0000000..49c6f0d
--- /dev/null
+++ b/src/ImaginaryOptics.jl
@@ -0,0 +1,41 @@
+"""
+    ImaginaryOptics.jl
+
+The code uses the negative sign convention for the propagation direction. Thus the gain and losses are defined using
+Thus a value of ϵ″>0 indicates gain and a value of ϵ″<0 indicates losses. The code uses the following convention for the sign of the imaginary part of the dielectric function:
+```math
+ϵ=ϵ′+im*ϵ″
+```
+"""
+module ImaginaryOptics
+import Unitful: μ0,ϵ0,c0,μm,ns,μs,s,cm,V,rad
+import Unitful.Length as UnitLength
+import Unitful.Time as UnitTime
+import Unitful: ustrip,unit, @u_str,𝐋,dimension
+using SparseArrays
+import Arpack: eigs # for eigenvalues
+@enum SOURCE_TYPE TFSF SOFT_SOURCE
+@enum MODE_TYPE TE TM
+@enum PML_TYPE PML_CLASSIC PML_RUMPF 
+@enum BOUNDARY_CONDITIONS begin
+    DIRICHLET
+    NEUMANN
+    ROBIN
+    PERIODIC
+    TRANSPARENT
+end
+
+@enum SOLVER_TYPE BEAM_PROPAGATION_METHODS DOMAIN_DECOMPOSITION_METHODS FDFD_METHODS
+abstract type AbstractField{T<:AbstractFloat} end
+abstract type AbstractSource{T<:AbstractFloat} end
+abstract type AbstractSolver{T<:AbstractFloat} end
+abstract type AbstractMaterial{T<:AbstractFloat} end
+include("utility.jl")
+include("common_solver_functions.jl")
+
+
+export μm,μ0,ϵ0,c0
+export propagation_phase
+# Write your package code here.
+
+end
\ No newline at end of file
diff --git a/src/beam_propagation_method.jl b/src/beam_propagation_method.jl
new file mode 100644
index 0000000..aec9511
--- /dev/null
+++ b/src/beam_propagation_method.jl
@@ -0,0 +1,52 @@
+
+
+
+""" 
+    BPM_solver{T}(Nx::Int,Ny::Int,dx_norm::T,dy_norm::T,ERzz::Matrix{Complex{T}};nonlinear::Bool=true,Npmly=10,kwargs...) where T<:AbstractFloat
+    BPM_solver{T}(NS::NTuple{2,Int},RES::NTuple{2,T},ERzz;kwargs...) where T<:AbstractFloat
+    BPM_solver{T}(NS::NTuple{2,Int},dx_norm::T,ERzz;kwargs...)       where T<:AbstractFloat
+    BPM_solver{T}(Nx::Int,Ny::Int,dx_norm::T,ERzz;kwargs...)         where T<:AbstractFloat
+
+    BPM_solver is a subtype of AbstractDomain that contains the information to solve the numerical problem. As the name says it is used to solve the problem using the beam
+    propagation method (BPM). The current version uses the nonlinear Crank Nicolson method to the direct problem and autodifferentiation to solve the inverse problem. The current
+    version only supports the paraxial approximation, future version will support Pade approximants to solve large angles.
+
+    # INPUTS
+    - `Nx::Int`: number of points in the x direction
+    - `Ny::Int`: number of points in the y direction
+    - `dx_norm::T`: grid spacing in the x direction in normalized units (i.e. dx_norm=dx*k₀)
+    - `dy_norm::T`: grid spacing in the y direction in normalized units (i.e. dy_norm=dy*k₀)
+    - `ERzz ::Matrix{Complex{T}}`: matrix containing the real part of the dielectric function. It is a misnomer, since in the case of active devices it 
+    can be complex and equal to ϵ_core-k_min^2.
+    - `nonlinear::Bool=false`: if true, the nonlinear term is included in the BPM equation (BPM has been made to solve the nonlinear optimization in 1D and 2D, by default it is true)
+    - `wpmly=10::Int=10`: number of PML points in the direction perpendicular to the propagation direction (y direction)
+    See also: [DomainFDFD](@ref), [DomainDecomposition](@ref), [DomainDecompositionLowMemory](@ref),
+"""
+struct BPM_solver{T} <: AbstractDomain{T}
+    nonlinear     :: Bool
+    N             :: Tuple{Int,Int}
+    laplacian     :: SparseMatrixCSC{Complex{T}}
+    epsilon       :: Matrix{T}
+    
+    function BPM_solver(Nx::Int,Ny::Int,dx_norm::T,dy_norm::T,epsilon::Vector{T};
+        nonlinear::Bool=true,wpmly=10,kwargs...) where T<:AbstractFloat   
+        
+        
+        @info("Nx: $Nx, Ny: $Ny, dx_norm: $dx_norm, dy_norm: $dy_norm, Npmly: $Npmly")
+        _,∇=Laplacian(1,Ny,dx_norm,dy_norm;wpmlx=0,wpmly=wpmly,check=2,kwargs...)
+        @debug("Nx: $Nx, Ny: $Ny, dx_norm: $dx_norm, dy_norm: $dy_norm, Npmly: $Npmly")
+        _,∇=Laplacian(1,Ny,dx_norm,dy_norm;wpmlx=0,wpmly=wpmly,check=2,kwargs...)
+        # Assume μ=1 that is the case for most materials
+        laplacian=∇[4]*∇[2]
+       
+        @debug("lapacialian created")
+        N=(Nx,Ny)
+        epsilon=reshape(epsilon,Nx,Ny)
+        new{T}(nonlinear,(Nx,Ny),laplacian,epsilon)
+    end
+end
+
+BPM_solver(NS::NTuple{2,Int},RES::NTuple{2,T},epsilon_real::Union{VecOrMat{T},VecOrMat{Complex{T}}};kwargs...) where T<:AbstractFloat =BPM_solver(NS...,RES...,epsilon_real;kwargs...)
+BPM_solver(NS::NTuple{2,Int},dx_norm::T,epsilon_real::Union{VecOrMat{T},VecOrMat{Complex{T}}};kwargs...)       where T<:AbstractFloat =BPM_solver(NS...,dx_norm,dx_norm,epsilon_real;kwargs...)
+BPM_solver(Nx::Int,Ny::Int,dx_norm::T,epsilon_real::Union{VecOrMat{T},VecOrMat{Complex{T}}};kwargs...)         where T<:AbstractFloat =BPM_solver(Nx,Ny,dx_norm,dx_norm,epsilon_real;kwargs...)
+
diff --git a/src/common_solver_functions.jl b/src/common_solver_functions.jl
new file mode 100644
index 0000000..ccfeb8d
--- /dev/null
+++ b/src/common_solver_functions.jl
@@ -0,0 +1,623 @@
+function yeeder2d_grid(ER2,UR2)
+    #Rumpf pp.132
+    # To remember, look at the matrix sx_ez
+    # E has 1,1 -> _ey has 1,2 -> _ex has 2,1 : the are the opposite of the z one in their comonent
+    # H has 2,2 -> _ey has 2,1 -> _ex has 1,2 : the are the opposite of the z one in their comonent
+   
+    Nx2,Ny2=size(ER2);
+    ERxx =ER2[2:2:Nx2,1:2:Ny2];
+    ERyy =ER2[1:2:Nx2,2:2:Ny2];
+    ERzz =ER2[1:2:Nx2,1:2:Ny2];
+    URxx =UR2[1:2:Nx2,2:2:Ny2];
+    URyy =UR2[2:2:Nx2,1:2:Ny2];
+    URzz =UR2[2:2:Nx2,1:2:Ny2];
+    return ERxx,ERyy,ERzz,URxx,URyy,URzz
+end
+
+"""
+yeeder2d(NS,RES_unit;periodic=false,kinc,k0);
+yeeder2d(NS,RES_norm;periodic=false,kinc);
+
+Creates the derivative matrices for the finite difference frequency domain method for the staggered Yee grid.
+    
+ # INPUT
+ 
+    - `NS::Tuple{Int,Int}`: number of points in the x and y direction
+    - `RES::Tuple{T,T}`: grid spacing in the x and y direction in normalized units (i.e. dx_norm=dx*k₀)
+
+=================
+NS    (Nx, Ny) Grid Size
+RES   (dx, dy) Grid Resolution
+BC    [xbc ybc] Boundary Conditions
+        0: Dirichlet boundary conditions
+        1: Periodic boundary conditions
+kinc  [kx ky] Incident Wave Vector
+      This argument is only needed for PBCs.
+%
+Note: For normalized grids use k0*RES and kinc/k0
+%
+OUTPUT ARGUMENTS
+=================
+DEX   Derivative Matrix wrt x for Electric Fields
+DEY   Derivative Matrix wrt to y for Electric Fields
+DHX   Derivative Matrix wrt to x for Magnetic Fields
+DHY   Derivative Matrix wrt to y for Magnetic Fields
+"""
+function yeeder2d(NS::Tuple{Int,Int},RES_norm::Tuple{T,T};boundary_conditions::BOUNDARY_CONDITIONS=DIRICHLET) where T<:AbstractFloat
+    @assert(NS[1]>=1 && NS[2]>=1, "NS must be greater or equal to 1")
+    @debug("use of yeeder2d unitless")
+    RES2_norm=(RES_norm[1]/2,RES_norm[2]/2)
+    # Build DEX
+    M=prod(NS);
+    # center diagonal
+    d0=-ones(M);
+    d1= ones(M-1);
+    d1[NS[1]:NS[1]:M-1].=0.0;
+
+    
+    (NS[1]>1) ? DEX=spdiagm(0=>d0,1=>d1) : DEX=spzeros(M,M);
+    d1= ones(M-NS[1]);
+    
+    (NS[2]>1) ? DEY=spdiagm(0=>d0,NS[1]=>d1) : DEY=spzeros(M,M);
+
+    return (DEX./RES_norm[1],DEY./RES_norm[2], -1 .*transpose(copy(DEX))./RES_norm[1], -1 .*transpose(copy(DEY))./RES_norm[2]) 
+end
+
+function yeeder2d(NS::Tuple{Int,Int},RES_unit::Tuple{UnitLength,UnitLength};λ₀=1.55μm,theta=0)
+    @debug("use of yeeder2d with units")
+    λ₀=oftype(RES_unit[1],λ₀)
+    k₀=2*π/λ₀;          # 1/𝐋
+    
+    RES_norm=(k₀*RES_unit[1],k₀*RES_unit[2])
+
+    return yeeder2d(NS,RES_norm)
+
+    # Build DEX
+    M=prod(NS);
+    # center diagonal
+    d0=-ones(M);
+    d1= ones(M-1);
+    d1[NS[1]:NS[1]:M-1].=0.0;
+
+    (NS[1]>1) ? DEX=spdiagm(0=>d0,1=>d1) : DEX=spzeros(M,M);
+
+    # Build DEY
+    d1= ones(M-NS[1]);
+    
+    (NS[2]>1) ? DEY=spdiagm(0=>d0,NS[1]=>d1) : DEY=spzeros(M,M);
+
+    return (DEX./RES[1],DEY./RES[2], -1 .*transpose(copy(DEX))./RES[1], -1 .*transpose(copy(DEY))./RES[2]) 
+end
+yeeder2d(Nx,Ny,dx,dy;λ₀=1.55μm)=yeeder2d((Nx,Ny),(dx,dy);λ₀=1.55μm)
+yeeder2d(Nx,Ny,dx,dy)=yeeder2d((Nx,Ny),(dx,dy))
+
+"""
+    Laplacian(Nx::Int,Ny::Int,dx::T,dy::T;μ_xx=nothing,μ_yy=nothing,ϵ_xx=nothing,ϵ_yy=nothing,wpmlx=0,wpmly=0,mode_type::MODE_TYPE=TE,use_pml=true,check=true, kwargs...) where T<:Union{AbstractFloat,Unitful.Length}
+
+    Laplacian is a function that creates the laplacain operators for the finite difference frequency domain method. It is a wrapper for the yeeder2d function that creates the curl operators and the apply_scpml_TE function that applies the perfectly matched layer to the curl operators. It also creates the laplacian operator for the finite difference frequency domain method.
+    The implementation is based on Rumpft work.
+    
+    # INPUTS
+    - `Nx::Int`: number of points in the x direction
+    - `Ny::Int`: number of points in the y direction
+    - `dx::T`: grid spacing in the x direction in normalized units (i.e. dx_norm=dx*k₀)
+    - `dy::T`: grid spacing in the y direction in normalized units (i.e. dy_norm=dy*k₀)
+    - `μ_xx::Vector{T}`: vector containing the real part of the magnetic permeability in the x direction. It is a misnomer, since in the case of active devices it
+    can be complex and equal to μ_core-k_min^2.
+    - `μ_yy::Vector{T}`: vector containing the real part of the magnetic permeability in the y direction.
+    - `ϵ_xx::Vector{T}`: vector containing the real part of the dielectric function in the x direction. 
+    - `ϵ_yy::Vector{T}`: vector containing the real part of the dielectric function in the y direction. It is a misnomer, since in the case of active devices it
+    can be complex and equal to ϵ_core-k_min^2.
+    - `wpmlx::Int=0`: number of PML points in the direction perpendicular to the propagation direction (x direction)
+    - `wpmly::Int=0`: number of PML points in the direction perpendicular to the propagation direction (y direction)
+    - `mode_type::MODE_TYPE=TE`: type of mode to be solved, TE or TM
+    - `use_pml::Bool=true`: if true, the perfectly matched layer is applied to the curl operators
+    - `check::Bool=true`: if true, the function checks if the number of PML points is zero, if it is zero it will throw an error, if false it will ignore the error
+
+    # OUTPUTS
+    - `∇²::SparseMatrixCSC{Complex{T}}`: laplacian operator
+    - `∇ᵉ_x::SparseMatrixCSC{Complex{T}}`: gradient operator in the x direction for the electric field
+    - `∇ᵉ_y::SparseMatrixCSC{Complex{T}}`: gradient operator in the y direction for the electric field
+    - `∇ʰ_x::SparseMatrixCSC{Complex{T}}`: gradient operator in the x direction for the magnetic field
+    - `∇ʰ_y::SparseMatrixCSC{Complex{T}}`: gradient operator in the y direction for the magnetic field
+
+
+"""
+function Laplacian(Nx,Ny,dx::T,dy::T;μ_xx=nothing,μ_yy=nothing,ϵ_xx=nothing,ϵ_yy=nothing,wpmlx=0,wpmly=0,mode_type::MODE_TYPE=TE,use_pml=true,check=true, kwargs...) where T<:Union{AbstractFloat}
+    if (wpmlx==0 || wpmly==0) 
+        
+        if (check==true) 
+            @error "Npmlx or Npmly is zero, this will result in a non-physical PML
+            set check=false to ignore this error" 
+            return nothing
+        else
+            @warn "Npmlx or Npmly is zero, this will result in a non-physical PML
+            Npmlx: $(wpmlx) Npmly: $(wpmly)"
+        end
+    end
+    NS=(Nx,Ny)
+    RES=(dx,dy) 
+    @debug ("NS: $(NS), RES: $(RES)")
+    @debug("NPML: $((wpmlx,wpmlx,wpmly,wpmly))")
+    ∇ᵉ_x,∇ᵉ_y,∇ʰ_x,∇ʰ_y=yeeder2d(Nx,Ny,dx,dy)
+    if use_pml
+        (∇ᵉ_x,∇ᵉ_y,∇ʰ_x,∇ʰ_y)=apply_scpml_TE(∇ᵉ_x,∇ᵉ_y,∇ʰ_x,∇ʰ_y,Nx,Ny,(dx,dy);NPML=(wpmlx,wpmlx,wpmly,wpmly),kwargs...);
+    end
+
+
+    if (isnothing(μ_xx) && isnothing(μ_yy) && isnothing(ϵ_xx) && isnothing(ϵ_yy))
+        @warn "No material properties specified, the laplacian will not be created only the curl operators"
+        ∇²=nothing
+        return (∇²,∇ᵉ_x,∇ᵉ_y,∇ʰ_x,∇ʰ_y)
+    end
+       
+    if (mode_type==TE || mode_type==H) 
+        @assert μ_xx!=nothing && μ_yy!=nothing "μ_xx and μ_yy must be specified for TE mode"
+        ONE=one(typeof(μ_yy[1]))
+        @debug "∇²: $(size(∇ᵉ_x)) $(size(∇ʰ_x)) $(size(∇ᵉ_y)) $(size(∇ʰ_y)) μ_xx: $(size(μ_xx)) μ_yy: $(size(μ_yy))"
+        ∇²=∇ʰ_x*spdiagm(0=>ONE./μ_yy[:])*∇ᵉ_x.+∇ʰ_y*spdiagm(0=>ONE./μ_xx[:])*∇ᵉ_y
+        return (∇²,(∇ᵉ_x,∇ᵉ_y,∇ʰ_x,∇ʰ_y))
+    elseif (mode_type==TM || mode_type==E)
+        @assert ϵ_xx!=nothing && ϵ_yy!=nothing "ϵ_xx and ϵ_yy must be specified for TM mode"
+        ONE=one(typeof(ϵ_xx[1]))
+        ∇²=∇ᵉ_x*spdiagm(0=>ONE./ϵ_yy[:])*∇ʰ_x.+DEY*spdiagm(0=>ONE./ϵ_xx[:])*∇ʰ_y
+        return ∇²,∇ᵉ_x,∇ᵉ_y,∇ʰ_x,∇ʰ_y
+    end
+
+    error("Mode has to be TE(H), TM(E)")
+end
+
+
+"""
+    SCPML_classic(a_max,σ_max,p,i,Nx)
+Classic implementation of stretched coordinates perfectly matched layer (SCPML) for the finite difference frequency domain method (FDFD).
+```math
+
+σ=σ_max(i/Nx)^p
+
+s=a_max(1+im*σ)
+```
+# Arguments
+- `a_max::Float`: maximum value of the stretching factor
+- `σ_max::Float`: maximum value of the conductivity
+- `p::Float`: stretching factor
+- `i::Int`: current position in the grid
+"""
+function SCPML_classic(a_max,σ_max,p,i,Nx)
+    σ=σ_max*(i/Nx)^p;
+    return a_max*(1-im*σ);
+end
+
+"""
+    SCPML_rumpf(a_max,σ_max,p,i,Nx)
+Implementation of stretched coordinates perfectly matched layer (SCPML) for the finite difference frequency domain method (FDFD) based on the work of Rumpft
+```math
+
+a=1+(a_max-1)*(i/Nx)^p
+σ=σ_max*sin(π*i/(2*Nx))^2
+
+s=a_max(1+im*σ*60*σ)
+```
+# Arguments
+- `a_max::Float`: maximum value of the stretching factor
+- `σ_max::Float`: maximum value of the conductivity
+- `p::Float`: stretching factor
+- `i::Int`: current position in the grid
+"""
+function SCPML_rumpf(a_max,σ_max,p,i,Nx)
+    a=1+(a_max-1)*(i/Nx)^p;
+    σ=σ_max*sin(π*i/(2*Nx))^2;
+
+    return a*(1-im*60*σ);    
+end
+
+function SC_PML(Nx::Int,Ny::Int,RES_norm::NTuple{2};
+    NPML::NTuple{4,Int}=(20,20,20,20),m=2.5,R=1e-7,λₘᵢₙ=1.55μm,
+    pml_type::PML_TYPE=PML_CLASSIC,kwargs...)    
+ 
+    sx2 = ones(Complex{Float32}, 2*Nx,2*Ny);
+    sy2 = ones(Complex{Float32}, 2*Nx,2*Ny);
+    (Nx2,Ny2)=2 .*(Nx,Ny);
+    NPML2=[2 .*npml for npml in NPML]
+    RES2_norm=(RES_norm[1]/2,RES_norm[2]/2)
+    ###########################################
+    d_PML=(NPML2[1]*RES2_norm[1],NPML2[2]*RES2_norm[1],NPML2[3]*RES2_norm[2],NPML2[4]*RES2_norm[2]);    
+    if pml_type==PML_CLASSIC
+        @debug("classic PML")
+        σₘₐₓ=((-(m+1)*log(R)./(2 .*d_PML))); # unitless
+        @info("σₘₐₓ: $(σₘₐₓ)")
+        aₘₐₓ=1.0
+        
+    elseif pml_type==PML_RUMPF
+        @debug("rumpf PML")
+        σₘₐₓ=[1.0,1.0,1.0,1.0]
+        aₘₐₓ=4.0
+    else
+        error("incorrect PML pml_type")
+    end
+    #return pml,σₘₐₓ
+    @debug("σₘₐₓ: $(σₘₐₓ)")
+    pml(σₘₐₓ,i,Nx)=(pml_type==PML_CLASSIC) ? SCPML_classic(aₘₐₓ,σₘₐₓ,m,i,Nx) :  SCPML_rumpf(aₘₐₓ,σₘₐₓ,m,i,Nx) 
+    
+    for nx=1:NPML2[1]
+        sx2[NPML2[1]-nx+1,:].=pml(σₘₐₓ[1],float(nx),float(NPML2[1]));
+    end
+    
+    for nx=1:NPML2[2]
+        sx2[Nx2-NPML2[2]+nx,:].=pml(σₘₐₓ[2],float(nx),float(NPML2[2]));
+    end
+    
+    for ny=1:NPML2[3]
+        sy2[:,NPML2[3]-ny+1].=pml(σₘₐₓ[3],float(ny),float(NPML2[3]));
+    end
+    
+    for ny=1:NPML2[4]
+        sy2[:,Ny2-NPML2[4]+ny,:].=pml(σₘₐₓ[4],float(ny),float(NPML2[4]));
+    end
+   # To remember, look at the matrix sx_ez
+   # E has 1,1 -> _ey has 1,2 -> _ex has 2,1 : the are the opposite of the z one in their comonent
+   # H has 2,2 -> _ey has 2,1 -> _ex has 1,2 : the are the opposite of the z one in their comonent
+   sx_ey = sx2[1:2:Nx2,2:2:Ny2]
+   sx_ez = sx2[1:2:Nx2,1:2:Ny2]
+   sy_ex = sy2[2:2:Nx2,1:2:Ny2]
+   sy_ez = sy2[1:2:Nx2,1:2:Ny2]
+   
+   sx_hy = sx2[2:2:Nx2,1:2:Ny2]
+   sx_hz = sx2[2:2:Nx2,2:2:Ny2]
+   sy_hx = sy2[1:2:Nx2,2:2:Ny2]
+   sy_hz = sy2[2:2:Nx2,2:2:Ny2]
+    #returns the stretching to be used to update both the matrix Dx and Dy AND the Robin BC
+    return  sx_ey,sx_ez,sy_ex,sy_ez,sx_hy,sx_hz,sy_hx,sy_hz #,sx,sy
+
+end
+
+
+
+function apply_scpml_TE(DEX,DEY,DHX,DHY,Nx::Int,Ny::Int,(dx,dy);kwargs...)
+    @debug "kwargs: $(keys(kwargs))"
+    (sx_ey,sx_ez,sy_ex,sy_ez,sx_hy,sx_hz,sy_hx,sy_hz)=SC_PML(Nx,Ny,(dx,dy);kwargs...)
+
+    DEX=spdiagm(0=>1.0 ./(sx_hy[:]))*DEX;
+    DEY=spdiagm(0=>1.0 ./(sy_hx[:]))*DEY;
+    DHX=spdiagm(0=>1.0 ./(sx_ez[:]))*DHX;
+    DHY=spdiagm(0=>1.0 ./(sy_ez[:]))*DHY;
+    return (DEX,DEY,DHX,DHY,dx,dy)
+end
+
+function apply_scpml_TM(DEX,DEY,DHX,DHY,Nx::Int,Ny::Int,(dx,dy);kwargs...)
+    
+    (sx_ey,sx_ez,sy_ex,sy_ez,sx_hy,sx_hz,sy_hx,sy_hz)=SC_PML(Nx,Ny,(dx,dy);kwargs...)
+
+    DEX=spdiagm(0=>1.0 ./(sx_hz[:]))*DEX;
+    DEY=spdiagm(0=>1.0 ./(sy_hz[:]))*DEY;
+    DHX=spdiagm(0=>1.0 ./(sx_ey[:]))*DHX;
+    DHY=spdiagm(0=>1.0 ./(sy_ex[:]))*DHY;
+    return (DEX,DEY,DHX,DHY,dx,dy)
+end
+
+#=
+function apply_scpml_TE(DEX,DEY,DHX,DHY,Nx::Int,Ny::Int,(dx,dy);kwargs...)
+    (_,sx_ez,_,sy_ez,sx_hy,_,sy_hx,_
+    )=SC_PML(Nx,Ny,(dx,dy);kwargs...)
+
+    DEX=spdiagm(0=>1.0 ./(sx_hy[:]))*DEX;
+    DEY=spdiagm(0=>1.0 ./(sy_hx[:]))*DEY;
+    DHX=spdiagm(0=>1.0 ./(sx_ez[:]))*DHX;
+    DHY=spdiagm(0=>1.0 ./(sy_ez[:]))*DHY;
+    return (DEX,DEY,DHX,DHY,dx,dy)
+end
+=#
+
+
+function slabmode(W::T,Ny,dy_norm::V; mode= 0, k₀=2π/1.55μm,n_core,n_cladding,n_substrate=nothing,shift=0.0μm) where {T,V<:AbstractFloat}
+    if isnothing(n_substrate);
+        n_substrate=n_cladding
+    end
+    
+    println("k₀   : $k₀");
+    NA=sqrt(n_core^2-n_cladding^2);  
+    @debug("k₀: $k₀ NA: $NA W: $W")
+    R=k₀*W*NA;  
+    M_allowed=Int(floor(2*R/pi));
+    println("Numerical Aperture   : $NA")
+    println("Normalized Frequency : $R")
+    println("Mode Allowed : $M_allowed")
+    N_used=Int(min(mode,M_allowed)+1);
+    
+    Ny_buffer=Int(floor(Ny/2));
+    Ny=Ny+2*Ny_buffer;
+    Ny2=Int(Ny*2);
+    dy2_norm=dy_norm/2;
+    dy2_unit=dy2_norm/k₀;
+    ER2=ones(1,Ny2).*n_cladding^2;
+    UR2=ones(1,Ny2);
+    Y2_wg=range(1,Ny2,Ny2).*dy2_unit;
+    Y2_wg= Y2_wg.-mean(Y2_wg);
+    
+    ER2= [abs(y.+shift).<= W/2 ? n_core.^2 : (((y.+shift).<-W/2) ? n_substrate.^2 : n_cladding.^2) for y in Y2_wg];
+    ER2=reshape(ER2,1,Ny2);
+
+    # EXTRACT YEE GRID ARRAYS
+    ERxx = ER2[1,1:2:Ny2]; ERxx=spdiagm(0=> ERxx[:]);
+    ERyy = ER2[1,2:2:Ny2]; ERyy=spdiagm(0=> ERyy[:]);
+    ERzz = ER2[1,1:2:Ny2]; ERzz=spdiagm(0=> ERzz[:]);
+    URxx = UR2[1,2:2:Ny2]; URxx=spdiagm(0=> URxx[:]);
+    URyy = UR2[1,1:2:Ny2]; URyy=spdiagm(0=> URyy[:]);
+    URzz = UR2[1,2:2:Ny2]; URzz=spdiagm(0=> URzz[:]);
+    
+
+    NS  = (1, Ny);
+    RES = (ones(V),dy_norm);
+    (DEX,DEY,DHX,DHY) = yeeder2d(NS,RES);
+    A = -(DHY/URzz*DEY + ERxx);
+    B = inv(Matrix(URyy));
+    ev      = -n_core.^2;
+    eigenvalues, eigenvectors = eigs(Matrix(A), B,which=:SR, check=2) #N_used+1,  which=:SR);
+    
+    D       = [sqrt(Complex{V}(x)) for x in eigenvalues];
+    NEFF    = -im*diagm(D);
+    
+    return real(NEFF)[:],eigenvectors[Ny_buffer+1:end-Ny_buffer,:],ERzz[Ny_buffer+1:end-Ny_buffer,:]
+
+end
+
+
+
+
+mutable struct Tensor
+    xx::Array;
+    yy::Array;
+    zz::Array;
+    xy::Array;
+    xz::Array;
+    yz::Array;
+end
+
+function create_FDFD2D(device,RES::NTuple{2};
+    n_core,n_cladding,λ=1.55μm,source_position,NPMLx=20,NPMLy=20,padding_x=0.0μm,padding_y=0.0μm)
+    #Nx=100;
+    #Ny=100;
+    #Δλ=0.01;
+    #λ0=1.55;
+    #device=Device(λ
+    dx=RES[1];
+    dy=RES[2];
+
+    NPML=(NPMLx,NPMLx,NPMLy,NPMLy);
+
+    xlim = (minimum(x for (_, (x, _)) in device.ports),maximum(x for (_, (x, _)) in  device.ports));
+    ylim = (minimum(y for ((_, y)) in device.polygon),maximum(y for ((_, y)) in  device.polygon));
+    
+    xlim=(xlim[1]-padding_x,xlim[2].+padding_x)
+    ylim=(ylim[1]-padding_y,ylim[2].+padding_y)
+    Nx=Int(Unitful.ustrip(ceil((xlim[2]-xlim[1])/dx)));
+    Ny=Int(Unitful.ustrip(ceil((ylim[2]-ylim[1])/dx)));
+
+    Nx2=2*Nx;
+    Ny2=2*Ny;
+    dx2=dx/2;
+    dy2=dy/2;
+
+    ER2=ones(Nx2,Ny2);
+    UR2=copy(ER2);
+
+    x_range=range(0.5,Nx-0.5)*dx;
+    y_range=range(0.5,Ny-0.5)*dy;
+    x_range=x_range.-mean(x_range);
+    y_range=y_range.-mean(y_range);
+
+    x2_range=range(0.5,Nx2-0.5)*dx2;
+    y2_range=range(0.5,Ny2-0.5)*dy2;
+    x2_range=x2_range.-mean(x2_range);
+    y2_range=y2_range.-mean(y2_range);
+
+    Y=[ (j) for  i in x_range,j in y_range];
+    X=[ (i) for i in x_range,j in y_range];
+
+
+    grid_points=[(x,y) for x in x2_range,y in y2_range];
+
+
+    # Define permittivities for inside and outside the polygon
+    epsilon_inside = n_core.^2;
+    epsilon_outside = n_cladding.^2;
+    pl=[ustrip.((x,y)) for (x,y) in device.polygon];
+
+# Assign permittivities based on whether the rotated points are inside or outside the polygon
+    for k in CartesianIndices(ER2)
+        (i,j)=Tuple(k)
+        point=grid_points[i,j];
+        ER2[i, j] = (inpolygon(ustrip.(point),pl)==1) ? epsilon_inside : epsilon_outside;
+    end
+
+    result=yeeder2d_grid(ER2,UR2)
+    ERzz=result[3];
+    URxx=result[4];
+    URyy=result[5];
+
+
+    ω=uconvert(u"rad/s",2*pi*c0/λ);
+    k₀=2π/λ;
+    println("create the FDFD matrix")
+    (DEX,DEY,DHX,DHY)=yeeder2d((Nx,Ny),(dx,dx).*k₀);
+    
+    println("Apply SCPML")
+    (DEX,DEY,DHX,DHY)=apply_scpml_TE(DEX,DEY,DHX,DHY,Nx,Ny,(dx,dx);NPML=NPML);
+    Dxx=DHX*spdiagm(0=>URyy[:])*DEX;
+    Dyy=DHY*spdiagm(0=>URxx[:])*DEY;
+
+    Ae=Dxx+Dyy+spdiagm(0=>real.(ERzz[:]));
+    println("Create source vector for TF/SF")
+    QQ=zeros(Nx,Ny);
+
+    QQ=X.<=source_position[1][1][1];
+
+    b0=spdiagm(0=>QQ[:])*Ae-Ae*spdiagm(0=>QQ[:]);
+
+    ER=Tensor(ERxx,ERyy,ERzz,[],[],[]);
+    UR=Tensor(URxx,URyy,URzz,[],[],[]);
+    RES=(dx*k₀,dy*k₀)
+    return Ae,b0, ER,UR,QQ, DEX,DEY,DHX,DHY,Nx,Ny,RES,X,Y,x_range,y_range
+end
+
+function create_optimization_elements(device)    
+    optimization_area=[];
+    if !isempty(device.monitor[OPTIMIZATION])
+        for (i,optim) in enumerate(device.monitor[OPTIMIZATION])
+            push!(optimization_area,
+            [(minimum(x for (x,_) in optim),(minimum(y for (_,y) in optim))),
+            (maximum(x for (x,_) in optim),(maximum(y for (_,y) in optim)))]);
+        end
+    end
+
+    source_position=[];
+    for source in device.monitor[SOURCE]
+        push!(source_position,[
+            (source[1][1],source[1][2]),(source[2][1],source[2][2])
+        ]);
+    end
+
+    measurment_position=[];
+    for measurment in device.monitor[MEASURMENT]
+        push!(measurment_position,[
+            (measurment[1][1],measurment[1][2]),(measurment[2][1],measurment[2][2])
+        ]);
+    end
+    return optimization_area,source_position,measurment_position    
+end
+
+"""
+    create_sources(W_wg_in,Nx,Ny,dy,position)
+
+    Create the sources for the TF/SF method and the modes for the optimization
+    # Arguments
+    * `W_wg_in::Float64`: width of the waveguide
+    * `X`: discretized x coordinates
+    * `Y`: discretized y coordinates
+    * `dx::Unitful.Length`: grid spacing in the y direction
+    * `dy::Unitful.Length`: grid spacing in the y direction
+    * `position`           : position of the source
+
+    # Returns
+    * `sources`            : source vector for TF/SF
+    * `mode`               : mode vector for optimization
+    * `phase`              : phase vector for optimization
+    * `eigenvalue`         : eigenvalue of the mode
+"""
+function create_sources(W_wg,X,Y,dx,dy,position;λ, noshift=false,kwargs...)
+    Nx= size(X,1)
+    Ny= size(X,2)
+    eigenvalue,e2,ycoor=slabmode(W_wg,Ny,dy;λ,kwargs...);
+    k₀=2π/λ;
+    sources = zeros(Complex{Float32},Nx, Ny, length(position));
+    phase   = zeros(Complex{Float32},Nx);
+    for (k, source)  in enumerate(position)
+        shifted_e2 = zeros(Complex{Float32}, size(e2));
+        if noshift
+            shifted_e2 = e2;
+        else
+            center = (source[2] .+ source[1]) ./ 2;
+            idx = findmin(abs.(Y[1, :] .- center[2]))[2];
+        
+            # Calculate the shift required for e2
+            local tmp_shift = idx - findmin(abs.(Y[1, :] .- 0μm))[2];
+            shifted_e2[max(1, 1+tmp_shift):min(end, end+tmp_shift)] = e2[max(1, 1-tmp_shift):min(end, end-tmp_shift)];
+        end
+        for i in axes(sources,1)
+            sources[i, :, k] = shifted_e2;
+        end
+    end
+    N=size(sources);#
+    modes=deepcopy(sources);
+    modes=reshape(modes,prod(N[1:2]),N[3]);
+
+    ddx=dx*k₀;
+    for i in axes(sources,1)
+        phase[i]=exp(-im*eigenvalue*(i-1)*ddx);
+    end
+    for k in CartesianIndices(sources)
+        i, j, l = Tuple(k)
+        sources[i, j, l] = sources[i, j, l] .* phase[i];
+    end
+
+    sources=reshape(sources,prod(N[1:2]),N[3]);
+
+    return sources,modes,phase,eigenvalue,center
+end
+
+function slabmode(W::T,Ny,dy_norm::V; mode= 0, k₀=2π/1.55μm,n_core,n_cladding,n_substrate=nothing,shift=0.0μm) where {T,V<:AbstractFloat}
+    if isnothing(n_substrate);
+        n_substrate=n_cladding
+    end
+    
+    println("k₀   : $k₀");
+    NA=sqrt(n_core^2-n_cladding^2);  
+    @debug("k₀: $k₀ NA: $NA W: $W")
+    R=k₀*W*NA;  
+    M_allowed=Int(floor(2*R/pi));
+    println("Numerical Aperture   : $NA")
+    println("Normalized Frequency : $R")
+    println("Mode Allowed : $M_allowed")
+    N_used=Int(min(mode,M_allowed)+1);
+    
+    Ny_buffer=Int(floor(Ny/2));
+    Ny=Ny+2*Ny_buffer;
+    Ny2=Int(Ny*2);
+    dy2_norm=dy_norm/2;
+    dy2_unit=dy2_norm/k₀;
+    #ER2=ones(1,Ny2).*n_cladding^2;
+    UR2=ones(Ny2);
+    Y2_wg=range(1,Ny2,Ny2).*dy2_unit;
+    Y2_wg= Y2_wg.-mean(Y2_wg);
+    
+    ER2= [abs(y.+shift).<= W/2 ? n_core.^2 : (((y.+shift).<-W/2) ? n_substrate.^2 : n_cladding.^2) for y in Y2_wg];
+    #ER2=reshape(ER2,1,Ny2);
+
+    # EXTRACT YEE GRID ARRAYS
+    ERxx = ER2[1:2:Ny2]; # ERxx=spdiagm(0=> ERxx[:]);
+    ERyy = ER2[2:2:Ny2]; #ERyy=spdiagm(0=> ERyy[:]);
+    ERzz = ER2[1:2:Ny2]; #ERzz=spdiagm(0=> ERzz[:]);
+    URxx = UR2[2:2:Ny2]; #URxx=spdiagm(0=> URxx[:]);
+    URyy = UR2[1:2:Ny2]; #URyy=spdiagm(0=> URyy[:]);
+    URzz = UR2[2:2:Ny2]; #URzz=spdiagm(0=> URzz[:]);
+    
+    NEFF,eigenvectors=slabmode(ERzz,URxx,URyy,buffer=0,dy_norm=dy_norm;mode=mode,n_core=n_core,n_cladding=n_cladding,n_substrate=n_substrate)
+    return NEFF,eigenvectors[Ny_buffer+1:end-Ny_buffer,:],ERzz[Ny_buffer+1:end-Ny_buffer,:]
+    
+end
+
+function slabmode(epsilon_zz::Vector{T},mu_xx::Vector{T},mu_yy::Vector{T},buffer::Int,dy_norm::T;
+    mode=0, n_core,n_cladding,n_substrate=nothing,k0=1.0) where T<:AbstractFloat
+    
+    ERzz=zeros(length(epsilon_zz)+2*buffer);
+    URxx=zeros(length(mu_xx)+2*buffer);
+    URyy=zeros(length(mu_yy)+2*buffer);
+    
+    @info("ERxx: $(size(ERyy)) URzz: $(size(URzz)) URyy: $(size(URxx))")
+
+    @info("epsilon_zz: $(size(epsilon_zz)) mu_yy: $(size(mu_yy)) mu_xx: $(size(mu_xx))")
+
+
+    ERzz[buffer+1:end-buffer]=epsilon_zz; ERzz[1:buffer].=ERzz[buffer+1]; ERzz[end-buffer+1:end].=ERzz[end-buffer];
+    
+    URyy[buffer+1:end-buffer]=mu_yy; URyy[1:buffer].=URyy[buffer+1]; URyy[end-buffer+1:end].=URyy[end-buffer];
+    URxx[buffer+1:end-buffer]=mu_xx; URxx[1:buffer].=URxx[buffer+1]; URxx[end-buffer+1:end].=URxx[end-buffer];
+
+    @info("ERxx: $(size(ERzz)) URzz: $(size(URyy)) URyy: $(size(URxx))")
+
+    Ny=length(ERzz);
+    NS  = (1, Ny);
+    RES = (oftype(dy_norm,1.0),dy_norm);
+    (DEX,DEY,DHX,DHY) = yeeder2d(NS,RES);
+    @info("DEX: $(size(DEX)) DEY: $(size(DEY)) DHX: $(size(DHX)) DHY: $(size(DHY)) URzz: $(size(URzz)) ERzz")
+    A = -(DHY/spdiagm(0=>URxx)*DEY + spdiagm(0=>ERzz.*k0^2));
+    B = inv(Matrix(spdiagm(0=>URyy)));
+    ev      = -n_core.^2;
+    eigenvalues, eigenvectors = eigs(Matrix(A), B,which=:SR, check=2) #N_used+1,  which=:SR);
+    
+    D       = [sqrt(Complex{T}(x)) for x in eigenvalues];
+    NEFF    = -im*diagm(D);
+    
+    return real(NEFF)[:],eigenvectors[buffer+1:end-buffer,mode+1] 
+
+end
\ No newline at end of file
diff --git a/src/domain_decomposition_method.jl b/src/domain_decomposition_method.jl
new file mode 100644
index 0000000..e69de29
diff --git a/src/finite_difference_frequency_domain.jl b/src/finite_difference_frequency_domain.jl
new file mode 100644
index 0000000..8b13789
--- /dev/null
+++ b/src/finite_difference_frequency_domain.jl
@@ -0,0 +1 @@
+
diff --git a/src/utility.jl b/src/utility.jl
new file mode 100644
index 0000000..dff7f55
--- /dev/null
+++ b/src/utility.jl
@@ -0,0 +1,27 @@
+
+function meshgrid(x_range,y_range)
+    Y=[ (j) for  i in x_range,j in y_range];
+    X=[ (i) for i in x_range,j in y_range];
+    return X,Y
+end      
+
+"""
+propagation_phase(n::T,dx_norm::T,dy_norm::T,Nx::Int,Ny::Int;theta=0) where T<:AbstractFloat
+propagation_phase(n::T,k0,dx_norm::M,dy_norm::M,Nx::Int,Ny::Int;;theta=0) where {T<:AbstractFloat,M}
+
+Computes the propagation phase for a given propagation constant n, the wavelength k0 and the length L. The phase is computed as exp(-im*n*L*cos(theta)) where theta is the angle between the propagation direction and the normal to the surface. If theta=0, the phase is exp(-im*n*L). If theta=pi/2, the phase is exp(-im*n*L*cos(pi/2))=exp(im*n*L*cos(pi/2))=exp(0)=1. If theta=pi, the phase is exp(-im*n*L*cos(pi))=exp(im*n*L*cos(pi))=exp(-im*n*L)=-1. If theta=3pi/2, the phase is exp(-im*n*L*cos(3pi/2))=exp(im*n*L*cos(3pi/2))=exp(0)=1. Thus the phase is exp(-im*n*L*cos(theta)).
+"""
+function propagation_phase(n::T,dx::T,dy::T,Nx::Int,Ny::Int;theta=0) where T<:AbstractFloat    
+    theta=oftype.(n,theta)
+    Lx=(0:dx:(Nx-1)*dx).*cos(theta)
+    Ly=0:dy:(Ny-1)*dy
+    Ly=reshape(Ly,1,Ny).*sin(theta)
+    return exp.(-im*n.*(Lx.+Ly))
+end
+function propagation_phase(n::T,k0,dx::M,dy::M,Nx::Int,Ny::Int;theta=0)  where {T<:AbstractFloat,M}
+    dx_norm=oftype.(n,k0.*dx)
+    dy_norm=oftype.(n,k0.*dy)
+    return propagation_phase(n,dx,dy,Nx,Ny;theta=theta)
+end
+
+
diff --git a/test/runtests.jl b/test/runtests.jl
new file mode 100644
index 0000000..446a3a5
--- /dev/null
+++ b/test/runtests.jl
@@ -0,0 +1,6 @@
+using ImaginaryOptics
+using Test
+
+@testset "ImaginaryOptics.jl" begin
+    # Write your tests here.
+end