Model of microscopic dissipation

Parameters

ω0 = 1;    
NBath = 300;   (* 
      Number of bath oscillators N *) 
ωmin = 0.1;  (* 
      Minimal frequency of bath oscillators *) 
ωmax = 2.1;  (* 
      Maximal frequency of bath oscillators *) 
ωi[
        i_] := ωmin + ((i - 1)/(NBath - 
                1)) (ωmax - ωmin);  (* 
      Frequency of the ith bath oscillator *) 
Coupling = 0.1;  (* 
      Coupling between the central oscillator and bath oscillators,  
      to be later divided by (Sqrt[NBath]) *) 
X0 = 
      1;  (* 
      Initial coordinate of the central oscillator *) 
tMax = 1500;  (* 
      Maximal time for plotting *) 
m = 
      1;  (* 
      Only this choice works for a moment in numerics since here M is a unit matrix *)

Solution

KMatr = 
      SparseArray[{{1, 1} :> ω(0^2), {i_, 
              i_} :> ωi[i - 1]^2, {i_, 1} :> 
            (Coupling/(Sqrt[NBath])), {1, i_} :> 
            (Coupling/(Sqrt[NBath]))}, {NBath + 1, 
          NBath + 1}]
MatrixForm[KMatr];   
Timing[ESK = Eigensystem[KMatr];]    
EValK = First[ESK];       
ΩΩ = (Sqrt[EValK]);   
EVecKHoriz = Last[ESK];    
AMatr = Transpose[EVecKHoriz];    
AMatrT = Transpose[AMatr]; (* This is a matrix of eigenvectors lying horizontally, our (A^(-1))=(A^T) *)
MatrixForm[AMatrT . AMatr];   

EE[
            t_] := ((((∑_(α = 1)^(NBath + 1))AMatr[[1, α]]^2 Cos[ΩΩ[[α]] t]))^2) + (1/ω)(0^2) ((((∑_(α = 1)^(NBath + 1))AMatr[[1, α]]^2 ΩΩ[[α]] Sin[ΩΩ[[α]] t]))^2)   

Γ0 = (π/4) ((((Coupling/(m ω0))))^2) 
        (1/(ωmax - ωmin));   (* 
      Analytical expression for damping *) 
EEAnalytical[t_] := 
      Exp[-2 Γ0 t];   (* 
      Analytical result for the normalized energy of the central oscillator *)

Plotting

PlotENum100 = 
      Plot[EE[t], {t, 0, tMax}, PlotRange -> {0, 1}, 
        PlotStyle -> Hue[0.7]]

PlotENum200 = 
      Plot[EE[t], {t, 0, tMax}, PlotRange -> {0, 1}]

PlotENum300 = 
      Plot[EE[t], {t, 0, tMax}, PlotRange -> {0, 1}, 
        PlotStyle -> Hue[0.3]]

PlotENum500 = 
      Plot[EE[t], {t, 0, tMax}, PlotRange -> {0, 1}, 
        PlotStyle -> Hue[0.4]]

PlotEAnal = 
      Plot[EEAnalytical[t], {t, 0, tMax}, PlotRange -> {0, 1}, 
        PlotStyle -> {Hue[0], AbsoluteThickness[2]}]

Show[PlotENum100, PlotENum200, PlotENum300, PlotEAnal]