Der geladene harmonische Oszillator

Picard-Lindelöf-Iteration

c0:	Anfangsvektor der Amplituden;
        N:	maximale Anzahl von Anregungen
        
        (Subscript[((c^(l + 1))), m]) = \(TraditionalForm\`\(c\^0\)\_m\ \ \  - \ 
        i/\(h\&~\)\ \ \(∫\_0\%t\ \ \ e\^\(i\ ω\_\(m\ n\)\ \ t'\)\ \ H\_\(m\ n\)\ \((t')\)\ dt'\)
\)
\(TraditionalForm\`ω\_\(m\ n\)\  = \ \((E\_m\ \  - \ 
              E\_n\ )\)/\(h\&~\)\)

\(\(Lindelöf[N_, ω1_, T_, tmax_, ymax_]\  := \ 
        Block[{ω0 = 1}, 
		c0\ \ \ \  = \ 
            Table[If[i == 0, 1, 0], {i, 0, 
                N}]; \ 
		c\ \ \ \ \ \  = \ 
            c0; 
Do[
		c\ \ \ \ \ \  = 
              Table[mi\  = \ 
                  m\  + \ 
                    1; \ \ c0[\([mi]\)] - \ I\ ω1 \((\ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ If[
                          m < N, \@\(m + 1\)\ \(∫\_0\%t\ \(E\^\(\(τ\)\(\ \)\(I\)\(\ \ \)\(ω0\)\(\ \)\)\) f[τ, T] \((c[\([mi + 
                                        1]\)] /. {t -> τ}\ )\)\ \ \ \[DifferentialD]τ\), 
                          0]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  + 
\ \ \ \ \ \ \ \ \ \ \ If[
                          m > 0, \(\@m\) \(∫\_0\%t\( E\^\(\(-τ\)\(\ \)\(I\)\(\ \)\(ω0\)\(\ \)\)\) f[τ, T] \((c[\([mi - 
                                        1]\)] /. {t -> τ}\ )\)\ \ \[DifferentialD]τ\), 
                          0]
	\ \ \ \ \ \ \ \ \ \ \ \ \ )\), {m, 0, N}]; 
		c\  = \ 
              Simplify[c]; 
Print[
              c[\([2]\)]]; 
		Plot[{Abs[c[\([1]\)]]\^2, 
                Abs[c[\([2]\)]]\^2, Abs[c[\([3]\)]]\^2}, {t, 0, tmax}, 
              PlotRange -> {{0, tmax}, {0, 
                    ymax}}, 
PlotStyle -> {{RGBColor[1, 0, 
                      0]}, {RGBColor[0, 1, 0]}, {RGBColor[0, 0, 
                      1]}}];
, {N}];
];\)\)

f[t_, 
          T_] := (E^(-t/T));
Lindelöf[5, ω1 =  .9, T =  .9, tmax = 10, 
        ymax = 1];

f[t_, T_] := 1;
Lindelöf[10, ω1 =  .5, T =  .9, tmax = 10, 
        ymax = 2];

f[t_, T_] := Cos[t/T];
Lindelöf[5, ω1 =  .2, T =  .9, tmax = 100, 
        ymax = 20];

f[t_] = ((((w1/2)))^2) ((((Cos[wp t]/wp  +  Cos[wm t]/wm  -  ((wm - wp)/(wm wp))))^(2 ))  +  (((Sin[wp t]/wp  +  Sin[wm t]/wm))^2))

w1  =  1;
wp  =  w0  +  w2;
wm  =  w0  -  w2;
w0  =  1;

g[t_] = f[t] /. w2 -> 1.1 // Simplify

Plot[f[t] /. w2 -> 0.9, {t, 0, 100}]