diff --git a/src/timecorrelations.jl b/src/timecorrelations.jl
index 80567240..99c2d90a 100644
--- a/src/timecorrelations.jl
+++ b/src/timecorrelations.jl
@@ -6,6 +6,8 @@ using QuantumOpticsBase
 using ..timeevolution, ..steadystate
 
 using FFTW
+using IterativeSolvers
+using LinearMaps
 
 
 """
@@ -157,5 +159,66 @@ function correlation2spectrum(tspan, corr; normalize_spec=false)
   omega, normalize_spec ? spec./maximum(spec) : spec
 end
 
+"""
+    timecorrelations.spectrum_iterative(omega, rho, op1, op2, H, J,
+                        method! = IterativeSolvers.bicgstabl!, args...;
+                        Jdagger=dagger.(J), rates=nothing, kwargs...)
+
+Iteratively solve the equation for the spectrum in Fourier space.
+"""
+function spectrum_iterative(omega, rho, op1, op2, H, J,
+                    method! = IterativeSolvers.bicgstabl!, args...;
+                    Jdagger=dagger.(J), rates=nothing, kwargs...)
+
+
+    bl = rho.basis_l
+    br = rho.basis_r
+    M = length(bl)
+
+    # Cache stuff
+    rho_cache = copy(rho)
+    drho = copy(rho)
+    Jrho_cache = copy(rho)
+
+    # Check reducibility
+    isreducible = timeevolution.check_master(rho,H,J,Jdagger,rates)
+    if isreducible
+        Hnh = timeevolution.nh_hamiltonian(H,J,Jdagger,rates)
+        Hnhdagger = dagger(Hnh)
+        dmaster_ = (drho,rho) -> timeevolution.dmaster_nh!(drho, Hnh, Hnhdagger, J, Jdagger, rates, rho, Jrho_cache)
+    else
+        dmaster_ = (drho,rho) -> timeevolution.dmaster_h!(drho, H, J, Jdagger, rates, rho, Jrho_cache)
+    end
+
+    i = 1
+    # Linear mapping
+    function f!(y,x)
+        # Reshape
+        drho.data .= @views reshape(y, M, M)
+        rho_cache.data .= @views reshape(x, M, M)
+        # Apply function
+        dmaster_(drho,rho_cache)
+        # Recast data
+        @views y .= reshape(drho.data, M^2)
+        @. y = im*omega[i]*x - y
+        return y
+    end
+
+    lm = LinearMaps.LinearMap{eltype(rho)}(f!,M^2;ismutating=true,issymmetric=false,ishermitian=false,isposdef=false)
+
+    R0 = op1*rho
+    x0 = @views reshape(R0.data, M^2)
+    b = reshape(copy(R0.data), M^2)
+
+
+    spec = similar(omega)
+    while i <= length(omega)
+        method!(x0,lm,b,args...;kwargs...)
+        spec[i] = 2*real(expect(op2,R0))
+        i += 1
+    end
+    return spec
+end
+
 
 end # module
diff --git a/test/test_timecorrelations.jl b/test/test_timecorrelations.jl
index 0075d971..b92ec179 100644
--- a/test/test_timecorrelations.jl
+++ b/test/test_timecorrelations.jl
@@ -64,4 +64,25 @@ omegaFFT, SFFT = timecorrelations.correlation2spectrum(tspan, exp_values)
 tspan[5] = tspan[4]
 @test_throws ArgumentError timecorrelations.correlation2spectrum(tspan, exp_values)
 
+# Test iterative spectrum
+bfock = FockBasis(10)
+η = 1.5
+κ = 1.0
+
+a = destroy(bfock)
+H = η*(a + a')
+J = [sqrt(2κ)*a]
+
+ρs = steadystate.iterative(H, J)
+@test isapprox(expect(a'*a,ρs), η^2 / κ^2; atol=1e-2)
+
+ω = range(-15, 15, length=200)
+S = timecorrelations.spectrum_iterative(ω, ρs, a, a', H, J)
+
+fwhm_idx = findmin(abs.(S .- 0.5*maximum(S)))[2]
+hwhm = abs(ω[fwhm_idx])
+@test abs(κ - hwhm) < 0.3
+
+# TODO: better tests and test lazy operators
+
 end # testset