diff --git a/ChainOhmT/mcdis2.jl b/ChainOhmT/mcdis2.jl
index af80714..4ce8834 100644
--- a/ChainOhmT/mcdis2.jl
+++ b/ChainOhmT/mcdis2.jl
@@ -69,7 +69,7 @@ MCDIS Multiple-component discretization procedure.
     in multi-component form.
 """
 
-function mcdis(N,eps0,quad::Function,Nmax)
+function mcdis(N,eps0,quad::Function,Nmax,idelta,mc,AB,wf,mp,irout)
 
     suc = true
     f = "Ncap exceeds Nmax in mcdis with irout = "
@@ -104,7 +104,7 @@ function mcdis(N,eps0,quad::Function,Nmax)
       xwm = Array{Float64}(undef, mc*Ncap, 2)
       for i=1:mc
         im1tn = (i-1)*Ncap
-        xw = quad(Ncap,i,uv)
+        xw = quad(Ncap,i,uv,mc,AB,wf)
         xwm[im1tn+1:im1tn+Ncap,:] = xw
       end
 
diff --git a/ChainOhmT/quadohmT.jl b/ChainOhmT/quadohmT.jl
index 2564827..9636388 100644
--- a/ChainOhmT/quadohmT.jl
+++ b/ChainOhmT/quadohmT.jl
@@ -1,73 +1,73 @@
-function quadfinT(N,i,uv)
+function quadfinT(N,i,uv,mc,AB,wf)
 
     if (i>1 && i<mc)
-        xw = trr(uv,i)
+        xw = trr(uv,i,AB,wf)
         return xw
     end
 
     if mc==1
       if (AB[i,1]!=-Inf)&&(AB[i,2]!=Inf)
-        xw = trr(uv,i)
+        xw = trr(uv,i,AB,wf)
         return xw
       elseif AB[i,1]!=-Inf
-          xw = rtr(uv,i)
+          xw = rtr(uv,i,AB,wf)
           return xw
       elseif AB[i,2]!=Inf
-          xw = ltr(uv,i)
+          xw = ltr(uv,i,AB,wf)
           return xw
       else
-          xw = str(uv,i)
+          xw = str(uv,i,wf)
           return xw
       end
     else
       if ((i==1 && AB[i,1]!=-Inf)||(i==mc && AB[i,2]!=Inf))
-        xw = trr(uv,i)
+        xw = trr(uv,i,AB,wf)
         return xw
       elseif i==1
-        xw = ltr(uv,i)
+        xw = ltr(uv,i,AB,wf)
         return xw
       end
     end
 
-    xw = rtr(uv,mc)
+    xw = rtr(uv,mc,AB,wf)
     return xw
 end
 
-function trr(t,i)
+function trr(t,i,AB,wf)
     s = Array{Float64}(undef, size(t))
     s[:,1] = ((AB[i,2]-AB[i,1])*t[:,1] .+ AB[i,2] .+ AB[i,1])./2
     s[:,2] = (AB[i,2]-AB[i,1]).*t[:,2].*wf(s[:,1],i)./2
     return s
 end
 
-function str(t,i)
+function str(t,i,wf)
     s = Array{Float64}(undef, size(t))
     s[:,1] = t[:,1]./(1-t[:,1].^2)
     s[:,2] = t[:,2].*wf(s[:,1],i).*(1+t[:,1].^2)./((1-t[:,1].^2).^2)
     return s
 end
 
-function rtr(t,i)
+function rtr(t,i,AB,wf)
     s = Array{Float64}(undef, size(t))
     s[:,1] = AB[i,1] .+ (t[:,1] .+ 1)./(-t[:,1] .+ 1)
     s[:,2] = 2*t[:,2].*wf(s[:,1],i)./((1 .-t[:,1]).^2)
     return s
 end
 
-function ltr(t,i)
+function ltr(t,i,AB,wf)
     s = Array{Float64}(undef, size(t))
     s[:,1] = -(-t[:,1].+1)./(t[:,1].+1) .+ AB[i,2]
     s[:,2] = 2*t[:,2].*wf(s[:,1],i)./((t[:,1] .+ 1).^2)
     return s
 end
 
-function wf(x,i)
+function ohmicspectraldensity_finiteT(x,i,α,s,ωc,β)
     if i==1
         y = 0
     elseif i==2
-        y = -( 2*a*abs.(x).^s ./ wc^(s-1)) .* ((coth.((beta/2).*x) .+ 1)./2) #.* exp(-abs(x)/wc)
+        y = -2.*( α*abs.(x).^s ./ ωc^(s-1)) .* (coth.((β/2).*x) .+ 1)./2 #.* exp(-abs(x)/wc)
     elseif i==3
-        y = ( 2*a*abs.(x).^s ./ wc^(s-1)) .* ((coth.((beta/2).*x) .+ 1)./2) #.* exp(-abs(x)/wc)
+        y = 2.*( α*abs.(x).^s ./ ωc^(s-1)) .* (coth.((β/2).*x) .+ 1)./2 #.* exp(-abs(x)/wc)
     elseif i==4
         y = 0
     end
diff --git a/ChainOhmT/stieltjes.jl b/ChainOhmT/stieltjes.jl
index d402884..8b2ad4e 100644
--- a/ChainOhmT/stieltjes.jl
+++ b/ChainOhmT/stieltjes.jl
@@ -9,8 +9,8 @@
     coefficients in the second column, of the nx2 array ab.
 """
 function stieltjes(N,xw) #return ab
-    tiny = 10*realmin
-    huge = .1*realmax
+    #tiny = 10*realmin
+    #huge = .1*realmax
 
     ## Remove data with zero weights ##
 
@@ -19,7 +19,7 @@ function stieltjes(N,xw) #return ab
     xw = xw[I,:]
     index = findfirst(x->x!=0, xw[:,2]) #index = min(find(xw(:,2)~=0))
     xw = xw[index:length(xw[:,2]),:]
-    Ibis = sortper(xw[:,1]) #xw = sortrows(xw,1)
+    Ibis = sortperm(xw[:,1]) #xw = sortrows(xw,1)
     xw = xw[Ibis,:]
     Ncap = size(xw,1)
 
@@ -27,7 +27,8 @@ function stieltjes(N,xw) #return ab
         error("N in sti out of range")
     end
 
-    s0 = ones(1,Ncap)*xw[:,2]
+    s0 = (ones(1,Ncap)*xw[:,2])[1]
+    ab = zeros(N, 2)
     ab[1,1] = transpose(xw[:,1])*xw[:,2]/s0
     ab[1,2] = s0
     if N==1
@@ -50,23 +51,22 @@ function stieltjes(N,xw) #return ab
     # % end
     # %
 
-    c = 1e240
-    p2 = c*p2
-    s0 = c^2*s0
+    #c = 1e240
+    #p2 = c*p2
+    #s0 = c^2*s0
 
     for k=1:N-1
       p0 = p1
       p1 = p2
-      p2 = (xw[:,1] - ab[k,1]).*p1 - ab[k,2]*p0
-      s1 = transpose(xw[:,2])*(p2.^2)
-      s2 = transpose(xw[:,1])*(xw[:,2].*(p2.^2))
-
-      if (max(abs(p2))>huge)||(abs(s2)>huge)
-        error("impending overflow in stieltjes for k=$(k)")
-      end
-      if abs(s1)<tiny
-        error("impending underflow in stieltjes for k=$(k)")
-      end
+      p2 = (xw[:,1] .- ab[k,1]).*p1 .- ab[k,2].*p0
+      s1 = (transpose(xw[:,2])*(p2.^2))[1]
+      s2 = (transpose(xw[:,1])*(xw[:,2].*(p2.^2)))[1]
+      # if (max(abs(p2))>huge)||(abs(s2)>huge)
+      #   error("impending overflow in stieltjes for k=$(k)")
+      # end
+      # if abs(s1)<tiny
+      #   error("impending underflow in stieltjes for k=$(k)")
+      # end
       ab[k+1,1] = s2/s1
       ab[k+1,2] = s1/s0
       s0 = s1
diff --git a/src/finitetemperature.jl b/src/finitetemperature.jl
new file mode 100644
index 0000000..1c797cb
--- /dev/null
+++ b/src/finitetemperature.jl
@@ -0,0 +1,114 @@
+using HDF5
+include("../ChainOhmT/quadohmT.jl")
+include("../ChainOhmT/mcdis2.jl")
+
+
+"""
+    chaincoeffs_finiteT(nummodes, β, ohmic=true; α, s, J, ωc=1, mc=4, mp=0, AB=nothing, iq=1, idelta=2, procedure=:Lanczos, Mmax=5000, save=true)
+
+Generate chain coefficients ``[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]`` for a harmonic bath at the inverse temperature β.
+
+By default a Ohmic spectral density ``J(ω) = 2αω_c (\\frac{ω}{ω_c})^s θ(ω-ω_c)`` is considered.
+Users can provide their own spectral density.
+
+# Arguments
+* nummodes: Number of bath modes
+* β: inverse temperature
+* ohmic: true if the spectral density is Ohmic, false if the user provides its own spectral density
+* α: Kondo parameter of the Ohmic spectral density 
+* s: ohmicity
+* J: user-provided spectral density. Should be a function f(x,i) where x is the frequency and i ∈ {1,...,mc} labels the intervals on which the SD is defined
+* ωc: the maximum frequency of the Ohmic spectral density
+* mc: the number of component intervals
+* mp: the number of points in the discrete part of the measure (mp=0 if there is none)
+* iq: a parameter to be set equal to 1, if the user provides his or her own quadrature routine, and different from 1 otherwise
+* idelta: a parameter whose default value is 1, but is preferably set equal to 2, if iq=1 and the user provides Gauss-type quadrature routines
+* procedure: choice between the Stieltjes and the Lanczos procedure
+* AB: component intervals
+* Mmax: maximum number of integration points
+* save: if true the coefficients are saved
+"""
+function chaincoeffs_finiteT(nummodes, β, ohmic=true; α=1, s=1, J=nothing, ωc=1, mc=4, mp=0, AB=nothing, iq=1, idelta=2, procedure=:Lanczos, Mmax=5000, save=true)
+
+    N = nummodes #Number of bath modes
+
+    if AB==nothing 
+        if mc==4
+            AB = [[-Inf -ωc];[-ωc 0];[0 ωc];[ωc Inf]]
+        else
+            throw(ArgumentError("An interval AB with mc = $mc components should have been provided."))
+        end
+    elseif length(AB) != mc
+        throw(ArgumentError("AB has a different number of intervals than mc = $mc."))             
+    end
+
+    if ohmic==true
+        wf(x,i) = ohmicspectraldensity_finiteT(x,i,α,s,ωc,β)
+    elseif J==nothing
+        throw(ArgumentError("A spectral density should have been provided."))
+    else
+        wf = J
+    end
+    
+    if procedure==:Lanczos  # choice between the Stieltjes (irout = 1) and the Lanczos procedure (irout != 1)
+        irout = 2 
+    elseif procedure==:Stieltjes
+        irout = 1 
+    else
+        throw(ArgumentError("Procedure should be either Lanczos or Stieltjes"))
+    end
+    
+    eps0=1e7*eps(Float64)
+
+    jacerg = zeros(N,2)
+
+    ab = 0.
+    ab, Mcap, kount, suc, uv = mcdis(N,eps0,quadfinT,Mmax,idelta,mc,AB,wf,mp,irout)
+    for m = 1:N-1
+        jacerg[m,1] = ab[m,1] #site energy e
+        jacerg[m,2] = sqrt(ab[m+1,2]) #hopping parameter t
+    end
+    jacerg[N,1] = ab[N,1]
+
+    eta = 0.
+    for i = 1:mc
+        xw = quadfinT(Mcap,i,uv,mc,AB,wf)
+        eta += sum(xw[:,2])
+    end
+    jacerg[N,2] = sqrt(eta) # system-chain coupling c
+
+    if save==true
+        # Write a HDF5 file
+        curdir = @__DIR__
+
+        if ohmic==true
+            Nstr=string(N)
+            astr=string(α)
+            sstr=string(s)
+            bstr=string(β)
+            # the "path" to the data inside of the h5 file is beta -> alpha -> s -> data (e, t or c)
+
+            # Write onsite energies
+            h5write("$curdir/ohmicT/chaincoeffs.h5", string("/N=",Nstr,"/α=",astr,"/s=",sstr,"/β=",bstr,"/e"), jacerg[1:N,1])
+            # Write hopping energies
+            h5write("$curdir/ohmicT/chaincoeffs.h5", string("/N=",Nstr,"/α=",astr,"/s=",sstr,"/β=",bstr,"/t"), jacerg[1:N-1,2])
+            # Write coupling coefficient
+            h5write("$curdir/ohmicT/chaincoeffs.h5", string("/N=",Nstr,"/α=",astr,"/s=",sstr,"/β=",bstr,"/c"), jacerg[N,2])
+
+        else
+            Nstr = string(N)
+            wstr = string(ωc)
+            bstr = string(β)
+            # the "path" to the data inside of the h5 file is N -> ωc -> beta -> data (e, t or c)
+
+            # Write onsite energies
+            h5write("$curdir/ohmicT/chaincoeffs.h5", string("/N=",Nstr,"/ωc=",wstr,"/β=",bstr,"/e"), jacerg[1:N,1])
+            # Write hopping energies
+            h5write("$curdir/ohmicT/chaincoeffs.h5", string("/N=",Nstr,"/ωc=",wstr,"/β=",bstr,"/t"), jacerg[1:N-1,2])
+            # Write coupling coefficient
+            h5write("$curdir/ohmicT/chaincoeffs.h5", string("/N=",Nstr,"/ωc=",wstr,"/β=",bstr,"/c"), jacerg[N,2])
+        end
+    end
+
+    return [jacerg[:,1], jacerg[1:N-1,2],jacerg[N,2]]
+end