function dcdt=bloch(t,c0)

N=(length(c0)-1)/2;
n=((-N):N)';
hbar = 1.024e-34;
lambda = 780e-9;
k=2*pi/lambda;
delta = 0*0.5*2*k;
g=9.81;



M = hbar*k^2/(2*2.36e4);
omegar = hbar*k^2/(2*M);
D=M*g/hbar*0.6e-3*0;
A = hbar/(2*M)*(2*n*k+delta+D).^2;
beta = 1*omegar;


%c0((2:length(c0))-1)

cn1 = [0 (c0((2:length(c0))-1))'];
cp1 = [(c0((1:(length(c0)-1))+1))' 0];

%whos

dcdt = -i*( A.*c0 + beta/2*(cn1' + cp1'));
%dcdt = dcdt/norm(dcdt);