diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml
new file mode 100644
index 00000000..604011f0
--- /dev/null
+++ b/.github/workflows/SpellCheck.yml
@@ -0,0 +1,13 @@
+name: Spell Check
+
+on: [pull_request]
+
+jobs:
+  typos-check:
+    name: Spell Check with Typos
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout Actions Repository
+        uses: actions/checkout@v4
+      - name: Check spelling
+        uses: crate-ci/typos@v1.27.0
\ No newline at end of file
diff --git a/.typos.toml b/.typos.toml
new file mode 100644
index 00000000..73648865
--- /dev/null
+++ b/.typos.toml
@@ -0,0 +1,2 @@
+[default.extend-words]
+ket = "ket"
\ No newline at end of file
diff --git a/README.md b/README.md
index f3bfceee..2661116b 100644
--- a/README.md
+++ b/README.md
@@ -65,7 +65,7 @@ QuantumToolbox.jl is equipped with a robust set of features:
 
 - **Quantum State and Operator Manipulation:** Easily handle quantum states and operators with a rich set of tools, with the same functionalities as QuTiP.
 - **Dynamical Evolution:** Advanced solvers for time evolution of quantum systems, thanks to the powerful [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl) package.
-- **GPU Computing:** Leverage GPU resources for high-performance computing. For example, you run the master equation direclty on the GPU with the same syntax as the CPU case.
+- **GPU Computing:** Leverage GPU resources for high-performance computing. For example, you run the master equation directly on the GPU with the same syntax as the CPU case.
 - **Distributed Computing:** Distribute the computation over multiple nodes (e.g., a cluster). For example, you can run hundreds of quantum trajectories in parallel on a cluster, with, again, the same syntax as the simple case.
 - **Easy Extension:** Easily extend the package, taking advantage of the Julia language features, like multiple dispatch and metaprogramming.
 
diff --git a/docs/src/index.md b/docs/src/index.md
index 4d1642c6..7447af11 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -14,7 +14,7 @@ QuantumToolbox.jl is equipped with a robust set of features:
 
 - **Quantum State and Operator Manipulation:** Easily handle quantum states and operators with a rich set of tools, with the same functionalities as QuTiP.
 - **Dynamical Evolution:** Advanced solvers for time evolution of quantum systems, thanks to the powerful [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl) package.
-- **GPU Computing:** Leverage GPU resources for high-performance computing. For example, you run the master equation direclty on the GPU with the same syntax as the CPU case.
+- **GPU Computing:** Leverage GPU resources for high-performance computing. For example, you run the master equation directly on the GPU with the same syntax as the CPU case.
 - **Distributed Computing:** Distribute the computation over multiple nodes (e.g., a cluster). For example, you can run undreds of quantum trajectories in parallel on a cluster, with, again, the same syntax as the simple case.
 - **Easy Extension:** Easily extend the package, taking advantage of the Julia language features, like multiple dispatch and metaprogramming.
 
diff --git a/docs/src/tutorials/lowrank.md b/docs/src/tutorials/lowrank.md
index 0a3c6c86..59e292bf 100644
--- a/docs/src/tutorials/lowrank.md
+++ b/docs/src/tutorials/lowrank.md
@@ -1,6 +1,6 @@
 # [Low rank master equation](@id doc-tutor:Low-rank-master-equation)
 
-In this tutorial, we will show how to solve the master equation using the low-rank method. For a detailed explaination of the method, we recommend to read the article [gravina2024adaptive](@cite).
+In this tutorial, we will show how to solve the master equation using the low-rank method. For a detailed explanation of the method, we recommend to read the article [gravina2024adaptive](@cite).
 
 As a test, we will consider the dissipative Ising model with a transverse field. The Hamiltonian is given by
 
diff --git a/docs/src/users_guide/states_and_operators.md b/docs/src/users_guide/states_and_operators.md
index 3cd71482..85993c2c 100644
--- a/docs/src/users_guide/states_and_operators.md
+++ b/docs/src/users_guide/states_and_operators.md
@@ -99,7 +99,7 @@ We can also create superpositions of states:
 ket = normalize(basis(5, 0) + basis(5, 1))
 ```
 
-where we have used the `normalize` function again to normalize the state. Apply the number opeartor again:
+where we have used the `normalize` function again to normalize the state. Apply the number operator again:
 
 ```@example states_and_operators
 n * ket
diff --git a/src/linear_maps.jl b/src/linear_maps.jl
index 8afd88e7..61ac4c96 100644
--- a/src/linear_maps.jl
+++ b/src/linear_maps.jl
@@ -18,7 +18,7 @@ A **linear map** is a transformation `L` that satisfies:
     L(cu) = cL(u)
     ```
 
-It is typically represented as a matrix with dimensions given by `size`, and this abtract type helps to define this map when the matrix is not explicitly available.
+It is typically represented as a matrix with dimensions given by `size`, and this abstract type helps to define this map when the matrix is not explicitly available.
 
 ## Methods
 
diff --git a/src/qobj/arithmetic_and_attributes.jl b/src/qobj/arithmetic_and_attributes.jl
index 96947bb2..5314d915 100644
--- a/src/qobj/arithmetic_and_attributes.jl
+++ b/src/qobj/arithmetic_and_attributes.jl
@@ -511,17 +511,17 @@ Quantum Object:   type=Operator   dims=[2]   size=(2, 2)   ishermitian=true
 function ptrace(QO::QuantumObject{<:AbstractArray,KetQuantumObject}, sel::Union{AbstractVector{Int},Tuple})
     _non_static_array_warning("sel", sel)
 
-    ns = length(sel)
-    if ns == 0 # return full trace for empty sel
+    n_s = length(sel)
+    if n_s == 0 # return full trace for empty sel
         return tr(ket2dm(QO))
     else
-        nd = length(QO.dims)
+        n_d = length(QO.dims)
 
-        (any(>(nd), sel) || any(<(1), sel)) && throw(
-            ArgumentError("Invalid indices in `sel`: $(sel), the given QuantumObject only have $(nd) sub-systems"),
+        (any(>(n_d), sel) || any(<(1), sel)) && throw(
+            ArgumentError("Invalid indices in `sel`: $(sel), the given QuantumObject only have $(n_d) sub-systems"),
         )
-        (ns != length(unique(sel))) && throw(ArgumentError("Duplicate selection indices in `sel`: $(sel)"))
-        (nd == 1) && return ket2dm(QO) # ptrace should always return Operator
+        (n_s != length(unique(sel))) && throw(ArgumentError("Duplicate selection indices in `sel`: $(sel)"))
+        (n_d == 1) && return ket2dm(QO) # ptrace should always return Operator
     end
 
     _sort_sel = sort(SVector{length(sel),Int}(sel))
@@ -534,17 +534,17 @@ ptrace(QO::QuantumObject{<:AbstractArray,BraQuantumObject}, sel::Union{AbstractV
 function ptrace(QO::QuantumObject{<:AbstractArray,OperatorQuantumObject}, sel::Union{AbstractVector{Int},Tuple})
     _non_static_array_warning("sel", sel)
 
-    ns = length(sel)
-    if ns == 0 # return full trace for empty sel
+    n_s = length(sel)
+    if n_s == 0 # return full trace for empty sel
         return tr(QO)
     else
-        nd = length(QO.dims)
+        n_d = length(QO.dims)
 
-        (any(>(nd), sel) || any(<(1), sel)) && throw(
-            ArgumentError("Invalid indices in `sel`: $(sel), the given QuantumObject only have $(nd) sub-systems"),
+        (any(>(n_d), sel) || any(<(1), sel)) && throw(
+            ArgumentError("Invalid indices in `sel`: $(sel), the given QuantumObject only have $(n_d) sub-systems"),
         )
-        (ns != length(unique(sel))) && throw(ArgumentError("Duplicate selection indices in `sel`: $(sel)"))
-        (nd == 1) && return QO
+        (n_s != length(unique(sel))) && throw(ArgumentError("Duplicate selection indices in `sel`: $(sel)"))
+        (n_d == 1) && return QO
     end
 
     _sort_sel = sort(SVector{length(sel),Int}(sel))
@@ -554,27 +554,27 @@ end
 ptrace(QO::QuantumObject, sel::Int) = ptrace(QO, SVector(sel))
 
 function _ptrace_ket(QO::AbstractArray, dims::Union{SVector,MVector}, sel)
-    nd = length(dims)
+    n_d = length(dims)
 
-    nd == 1 && return QO, dims
+    n_d == 1 && return QO, dims
 
-    qtrace = filter(i -> i ∉ sel, 1:nd)
+    qtrace = filter(i -> i ∉ sel, 1:n_d)
     dkeep = dims[sel]
     dtrace = dims[qtrace]
-    nt = length(dtrace)
+    n_t = length(dtrace)
 
     # Concatenate qtrace and sel without losing the length information
     # Tuple(qtrace..., sel...)
-    qtrace_sel = ntuple(Val(nd)) do i
-        if i <= nt
+    qtrace_sel = ntuple(Val(n_d)) do i
+        if i <= n_t
             @inbounds qtrace[i]
         else
-            @inbounds sel[i-nt]
+            @inbounds sel[i-n_t]
         end
     end
 
     vmat = reshape(QO, reverse(dims)...)
-    topermute = reverse(nd + 1 .- qtrace_sel)
+    topermute = reverse(n_d + 1 .- qtrace_sel)
     vmat = permutedims(vmat, topermute) # TODO: use PermutedDimsArray when Julia v1.11.0 is released
     vmat = reshape(vmat, prod(dkeep), prod(dtrace))
 
@@ -582,33 +582,33 @@ function _ptrace_ket(QO::AbstractArray, dims::Union{SVector,MVector}, sel)
 end
 
 function _ptrace_oper(QO::AbstractArray, dims::Union{SVector,MVector}, sel)
-    nd = length(dims)
+    n_d = length(dims)
 
-    nd == 1 && return QO, dims
+    n_d == 1 && return QO, dims
 
-    qtrace = filter(i -> i ∉ sel, 1:nd)
+    qtrace = filter(i -> i ∉ sel, 1:n_d)
     dkeep = dims[sel]
     dtrace = dims[qtrace]
-    nk = length(dkeep)
-    nt = length(dtrace)
-    _2_nt = 2 * nt
+    n_k = length(dkeep)
+    n_t = length(dtrace)
+    _2_n_t = 2 * n_t
 
     # Concatenate qtrace and sel without losing the length information
     # Tuple(qtrace..., sel...)
-    qtrace_sel = ntuple(Val(2 * nd)) do i
-        if i <= nt
+    qtrace_sel = ntuple(Val(2 * n_d)) do i
+        if i <= n_t
             @inbounds qtrace[i]
-        elseif i <= _2_nt
-            @inbounds qtrace[i-nt] + nd
-        elseif i <= _2_nt + nk
-            @inbounds sel[i-_2_nt]
+        elseif i <= _2_n_t
+            @inbounds qtrace[i-n_t] + n_d
+        elseif i <= _2_n_t + n_k
+            @inbounds sel[i-_2_n_t]
         else
-            @inbounds sel[i-_2_nt-nk] + nd
+            @inbounds sel[i-_2_n_t-n_k] + n_d
         end
     end
 
     ρmat = reshape(QO, reverse(vcat(dims, dims))...)
-    topermute = reverse(2 * nd + 1 .- qtrace_sel)
+    topermute = reverse(2 * n_d + 1 .- qtrace_sel)
     ρmat = permutedims(ρmat, topermute) # TODO: use PermutedDimsArray when Julia v1.11.0 is released
     ρmat = reshape(ρmat, prod(dkeep), prod(dkeep), prod(dtrace), prod(dtrace))
     res = map(tr, eachslice(ρmat, dims = (1, 2)))
diff --git a/src/qobj/boolean_functions.jl b/src/qobj/boolean_functions.jl
index 5262d7e9..869d1464 100644
--- a/src/qobj/boolean_functions.jl
+++ b/src/qobj/boolean_functions.jl
@@ -100,6 +100,6 @@ SciMLOperators.iscached(A::AbstractQuantumObject) = iscached(A.data)
 @doc raw"""
     SciMLOperators.isconstant(A::AbstractQuantumObject)
 
-Test whether the [`AbstractQuantumObject`](@ref) `A` is constant in time. For a [`QuantumObject`](@ref), this function returns `true`, while for a [`QuantumObjectEvolution`](@ref), this function returns `true` if the operator is contant in time.
+Test whether the [`AbstractQuantumObject`](@ref) `A` is constant in time. For a [`QuantumObject`](@ref), this function returns `true`, while for a [`QuantumObjectEvolution`](@ref), this function returns `true` if the operator is constant in time.
 """
 SciMLOperators.isconstant(A::AbstractQuantumObject) = isconstant(A.data)
diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index fbb51b8f..85491c6c 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -20,7 +20,7 @@ function _sprepost(A, B) # for any other input types
 end
 
 ## if input is AbstractSciMLOperator 
-## some of them are optimzed to speed things up
+## some of them are optimized to speed things up
 ## the rest of the SciMLOperators will just use lazy tensor (and prompt a warning)
 _spre(A::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_spre(A.A, Id))
 _spre(A::ScaledOperator, Id::AbstractMatrix) = ScaledOperator(A.λ, _spre(A.L, Id))
diff --git a/src/steadystate.jl b/src/steadystate.jl
index 097611d9..606b0907 100644
--- a/src/steadystate.jl
+++ b/src/steadystate.jl
@@ -106,11 +106,11 @@ function _steadystate(
     idx_range = collect(1:N)
     rows = _get_dense_similar(L_tmp, N)
     cols = _get_dense_similar(L_tmp, N)
-    datas = _get_dense_similar(L_tmp, N)
+    vals = _get_dense_similar(L_tmp, N)
     fill!(rows, 1)
     copyto!(cols, N .* (idx_range .- 1) .+ idx_range)
-    fill!(datas, weight)
-    Tn = sparse(rows, cols, datas, N^2, N^2)
+    fill!(vals, weight)
+    Tn = sparse(rows, cols, vals, N^2, N^2)
     L_tmp = L_tmp + Tn
 
     (haskey(kwargs, :Pl) || haskey(kwargs, :Pr)) && error("The use of preconditioners must be defined in the solver.")
@@ -160,11 +160,11 @@ function _steadystate(
     idx_range = collect(1:N)
     rows = _get_dense_similar(L_tmp, N)
     cols = _get_dense_similar(L_tmp, N)
-    datas = _get_dense_similar(L_tmp, N)
+    vals = _get_dense_similar(L_tmp, N)
     fill!(rows, 1)
     copyto!(cols, N .* (idx_range .- 1) .+ idx_range)
-    fill!(datas, weight)
-    Tn = sparse(rows, cols, datas, N^2, N^2)
+    fill!(vals, weight)
+    Tn = sparse(rows, cols, vals, N^2, N^2)
     L_tmp = L_tmp + Tn
 
     ρss_vec = L_tmp \ v0 # This is still not supported on GPU, yet
diff --git a/src/time_evolution/mcsolve.jl b/src/time_evolution/mcsolve.jl
index 31d5b1ec..07b66b51 100644
--- a/src/time_evolution/mcsolve.jl
+++ b/src/time_evolution/mcsolve.jl
@@ -38,14 +38,14 @@ function LindbladJumpAffect!(integrator)
     end
     cumsum!(cumsum_weights_mc, weights_mc)
     r = rand(traj_rng) * sum(weights_mc)
-    collaps_idx = getindex(1:length(weights_mc), findfirst(>(r), cumsum_weights_mc))
-    mul!(cache_mc, c_ops[collaps_idx], ψ)
+    collapse_idx = getindex(1:length(weights_mc), findfirst(>(r), cumsum_weights_mc))
+    mul!(cache_mc, c_ops[collapse_idx], ψ)
     normalize!(cache_mc)
     copyto!(integrator.u, cache_mc)
 
     random_n[] = rand(traj_rng)
     jump_times[internal_params.jump_times_which_idx[]] = integrator.t
-    jump_which[internal_params.jump_times_which_idx[]] = collaps_idx
+    jump_which[internal_params.jump_times_which_idx[]] = collapse_idx
     internal_params.jump_times_which_idx[] += 1
     if internal_params.jump_times_which_idx[] > length(jump_times)
         resize!(jump_times, length(jump_times) + internal_params.jump_times_which_init_size)
diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl
index a691ff0e..3d126b9f 100644
--- a/src/time_evolution/ssesolve.jl
+++ b/src/time_evolution/ssesolve.jl
@@ -5,11 +5,11 @@ export ssesolveProblem, ssesolveEnsembleProblem, ssesolve
 
 A struct to represent the diffusion operator. This is used to perform the diffusion process on N different Wiener processes.
 =#
-struct DiffusionOperator{T,OT<:Tuple{Vararg{AbstractSciMLOperator}}} <: AbstractSciMLOperator{T}
-    ops::OT
-    function DiffusionOperator(ops::OT) where {OT}
+struct DiffusionOperator{T,OpType<:Tuple{Vararg{AbstractSciMLOperator}}} <: AbstractSciMLOperator{T}
+    ops::OpType
+    function DiffusionOperator(ops::OpType) where {OpType}
         T = mapreduce(eltype, promote_type, ops)
-        return new{T,OT}(ops)
+        return new{T,OpType}(ops)
     end
 end
 
diff --git a/src/time_evolution/time_evolution_dynamical.jl b/src/time_evolution/time_evolution_dynamical.jl
index aa9660b0..270168df 100644
--- a/src/time_evolution/time_evolution_dynamical.jl
+++ b/src/time_evolution/time_evolution_dynamical.jl
@@ -8,12 +8,12 @@ function _reduce_dims(
     sel,
     reduce,
 ) where {T,N,DT<:Integer}
-    nd = length(dims)
+    n_d = length(dims)
     dims_new = zero(dims)
     dims_new[sel] .= reduce
     @. dims_new = dims - dims_new
 
-    if nd == 1
+    if n_d == 1
         ρmat = similar(QO, dims_new[1], dims_new[1])
         copyto!(ρmat, view(QO, 1:dims_new[1], 1:dims_new[1]))
     else
@@ -32,12 +32,12 @@ function _increase_dims(
     sel,
     increase,
 ) where {T,N,DT<:Integer}
-    nd = length(dims)
+    n_d = length(dims)
     dims_new = MVector(zero(dims)) # Mutable SVector
     dims_new[sel] .= increase
     @. dims_new = dims + dims_new
 
-    if nd == 1
+    if n_d == 1
         ρmat = similar(QO, dims_new[1], dims_new[1])
         fill!(selectdim(ρmat, 1, dims[1]+1:dims_new[1]), 0)
         fill!(selectdim(ρmat, 2, dims[1]+1:dims_new[1]), 0)
@@ -46,8 +46,8 @@ function _increase_dims(
         ρmat2 = similar(QO, reverse(vcat(dims_new, dims_new))...)
         ρmat = reshape(QO, reverse(vcat(dims, dims))...)
         for i in eachindex(sel)
-            fill!(selectdim(ρmat2, nd - sel[i] + 1, dims[sel[i]]+1:dims_new[sel[i]]), 0)
-            fill!(selectdim(ρmat2, 2 * nd - sel[i] + 1, dims[sel[i]]+1:dims_new[sel[i]]), 0)
+            fill!(selectdim(ρmat2, n_d - sel[i] + 1, dims[sel[i]]+1:dims_new[sel[i]]), 0)
+            fill!(selectdim(ρmat2, 2 * n_d - sel[i] + 1, dims[sel[i]]+1:dims_new[sel[i]]), 0)
         end
         copyto!(view(ρmat2, reverse!(repeat([1:n for n in dims], 2))...), ρmat)
         ρmat = reshape(ρmat2, prod(dims_new), prod(dims_new))