clear all;

color = 'b';

%------------------------------------------------------------

% Jaynes cummings oscillation with many oscillators

%------------------------------------------------------------

M = 100; % Number of harmonic oscillators



wa = 1; % atomic frequency

wk = linspace(0,4,M); % harmonic oscillator frequency

W = 0.1; % Rabi frequency, i.e. coupling coefficient



% Build the coupling hamiltonian

H = diag([wa wk]);

H(1,2:(1+M)) = W*ones(1,M);

H(2:(1+M),1) = W*ones(1,M);



% Initial state, |e,0,...,0>

psi_i = zeros(1+M,1);

psi_i(1) = 1;



% Solve the Hamiltonian

[V D] = eig(H);

eigbasis_i = V'*psi_i;



N = 100;

t = linspace(0,100,N);

prob = zeros(1,N);

for k = 1:N

    eigbasis_t = expm(-i*D*t(k))*eigbasis_i;

    psi_t = V*eigbasis_t;

    prob_t = abs(psi_t).^2;

    prob(k) = prob_t(1);

end

plot(t,prob,color); hold on;