NA = 300; NB = 100; % Nr of oscillators
qA = 500; qB = 1000; % Initial energy
q = qA + qB; % Total energy
N = NA + NB; % Total oscillators
state = zeros(N,1); % Microstate of the system
% state(1:NA) is part A, state(NA+1:NA+NB) is part B
% Generate initial, random microstate
placeA = randi(NA,qA,1);
figure(2)
for ip = 1:qA
    state(placeA(ip)) = state(placeA(ip)) + 1;
end
placeB = randi(NB,qB,1);
for ip = 1:qB
    i = NA + placeB(ip);
    state(i) = state(i) + 1;
end
% Simulate random state development
nstep = 40000; % Number of steps
EA = zeros(nstep,1); EB = zeros(nstep,1);
for istep = 1:nstep
    i1 = randi(N,1,1); % Select random osc. i1
    if (state(i1)>0) % Does it have any energy ?
        i2 = randi(N,1,1); % Select random osc. i2
        state(i2) = state(i2) + 1; % Transfer
        state(i1) = state(i1) - 1; % energy
    end
    EA(istep) = double(sum(state(1:NA)))/NA;
    EB(istep) = double(sum(state(NA+1:end)))/NB;
    if rem(istep,100)==1
        subplot(2,1,1) % Microstate of system
        plot((1:NA),state(1:NA),'b',(1:NB)+NA,state(NA+1:end),'r');
        xlabel('i'); ylabel('q_i'); drawnow
        subplot(2,1,2); % Avg energy in each system
        plot((1:istep),EA(1:istep),'-b',(1:istep),EB(1:istep),'-r');
        xlabel('t'); ylabel('q_A/N_A , q_B/N_B'), drawnow
    end
end
line([0 nstep],[q/N q/N])