diff --git a/examples/sbm_zero_temperature.jl b/examples/sbm_zero_temperature.jl
new file mode 100644
index 0000000..2218937
--- /dev/null
+++ b/examples/sbm_zero_temperature.jl
@@ -0,0 +1,83 @@
+"""
+    Example of a zero-temperature Spin-Boson Model with an hard cut-off Ohmic spectral density J(ω) = 2αω when ω < ωc and 0 otherwise
+
+    The dynamics is simulated using the T-TEDOPA method that maps the normal modes environment into a non-uniform tight-binding chain.
+
+    H = \frac{ω_0}{2} σ_z + Δ σ_x + c_0 σ_x(b_0^\dagger + b_0) + \sum_{i=0}^{N-1} t_i (b_{i+1}^\dagger b_i +h.c.) + \sum_{i=0}^{N-1} ϵ_i b_i^\dagger b_i  
+"""
+
+using MPSDynamics, Plots, LaTeXStrings
+
+#----------------------------
+# Physical parameters
+#----------------------------
+
+d = 6 # number of Fock states of the chain modes
+
+N = 30 # length of the chain
+
+α = 0.1 # coupling strength
+
+Δ = 0.0 # tunneling 
+
+ω0 = 0.2 # TLS gap
+
+s = 1 # ohmicity
+
+cpars = chaincoeffs_ohmic(N, α, s) # chain parameters, i.e. on-site energies ϵ_i, hopping energies t_i, and system-chain coupling c_0
+
+#-----------------------
+# Simulation parameters
+#-----------------------
+
+dt = 0.5 # time step
+
+tfinal = 30.0 # simulation time
+
+method = :TDVP1 # time-evolution method
+
+D = [2,4,6] # MPS bond dimension
+
+#---------------------------
+# MPO and initial state MPS
+#---------------------------
+
+H = spinbosonmpo(ω0, Δ, d, N, cpars) # MPO representation of the Hamiltonian
+
+ψ = unitcol(1,2) # Initial up-z system state 
+
+A = productstatemps(physdims(H), state=[ψ, fill(unitcol(1,d), N)...]) # MPS representation of |ψ>|Vacuum>
+
+#---------------------------
+# Definition of observables
+#---------------------------
+
+ob1 = OneSiteObservable("sz", sz, 1)
+
+ob2 = OneSiteObservable("chain mode occupation", numb(d), (2,N+1))
+
+ob3 = TwoSiteObservable("SXdisp", sx, disp(d), [1], collect(2:N+1))
+
+#-------------
+# Simulation
+#------------
+
+A, dat = runsim(dt, tfinal, A, H;
+                name = "ohmic spin boson model",
+                method = method,
+                obs = [ob2,ob3],
+                convobs = [ob1],
+                params = @LogParams(N, d, α, Δ, ω0, s),
+                convparams = D,
+                verbose = false,
+                save = true,
+                plot = true,
+                );
+
+#----------
+# Plots
+#----------
+
+plot(dat["data/times"], dat["convdata/sz"], label=["Dmax = 2" "Dmax = 4" "Dmax = 6"], xlabel=L"t",ylabel=L"\sigma_z")
+
+heatmap(dat["data/times"], collect(1:N), abs.(dat["data/SXdisp"][1,:,:]), xlabel=L"t",ylabel="chain mode")
diff --git a/src/models.jl b/src/models.jl
index 505e8af..6acdf72 100644
--- a/src/models.jl
+++ b/src/models.jl
@@ -7,7 +7,7 @@ dn(ops...) = permutedims(cat(reverse(ops)...; dims=3), [3,1,2])
 Return the MPO representation of the `N`-spin XYZ Hamiltonian with external field ``\\vec{h}=(h_x, 0, h_z)``.
 
 ``
-H = \\sum_{n=1}^{N-1} -J_x\\sigma_x^{n} \\sigma_x^{n+1} - J_y\\sigma_y^{n} \\sigma_y^{n+1} - J_z\\sigma_z^{n} \\sigma_z^{n+1} + \\sum_{n=1}^{N}(- h_x\\sigma_x^{n} - h_z\\sigma_z^{n})  
+H = \\sum_{n=1}^{N-1} -J_x σ_x^{n} σ_x^{n+1} - J_y σ_y^{n} σ_y^{n+1} - J_z σ_z^{n} σ_z^{n+1} + \\sum_{n=1}^{N}(- h_x σ_x^{n} - h_z σ_z^{n})  
 ``.
 
 """
@@ -345,7 +345,7 @@ end
 Generate MPO for a spin-1/2 coupled to a chain of harmonic oscillators, defined by the Hamiltonian
 
 ``
-H = \\frac{ω_0}{2}σ_z + Δσ_x + c_0σ_x(b_k^\\dagger+b_k) + \\sum_{i=0}^{N-1} t_i (b_{i+1}^\\dagger b_i +h.c.) + \\sum_{i=0}^{N-1} ϵ_ib_i^\\dagger b_i
+H = \\frac{ω_0}{2}σ_z + Δσ_x + c_0σ_x(b_0^\\dagger+b_0) + \\sum_{i=0}^{N-1} t_i (b_{i+1}^\\dagger b_i +h.c.) + \\sum_{i=0}^{N-1} ϵ_ib_i^\\dagger b_i
 ``.
 
 The spin is on site 1 of the MPS and the bath modes are to the right.