for j=1:3000
for i=1:200
[dbx1,dby1,dbz1]=nonlinear(bx,by,bz,kx,ky,kz,ksquare);
dbx1=dbx1-nu*ksquare.*bx;
dby1=dby1-nu*ksquare.*by;
dbz1=dbz1-nu*ksquare.*bz;
dbx1(2,1,1)=dbx1(2,1,1)-g/2;
dbx1(M,1,1)=dbx1(M,1,1)-g/2;
bx1=bx+0.5*deltat*dbx1;
by1=by+0.5*deltat*dby1;
bz1=bz+0.5*deltat*dbz1;

[dbx2,dby2,dbz2]=nonlinear(bx1,by1,bz1,kx,ky,kz,ksquare);
dbx2=dbx2-nu*ksquare.*bx1;
dby2=dby2-nu*ksquare.*by1;
dbz2=dbz2-nu*ksquare.*bz1;
dbx2(2,1,1)=dbx2(2,1,1)-g/2;
dbx2(M,1,1)=dbx2(M,1,1)-g/2;
bx2=bx+0.5*deltat*dbx2;
by2=by+0.5*deltat*dby2;
bz2=bz+0.5*deltat*dbz2;

[dbx3,dby3,dbz3]=nonlinear(bx2,by2,bz2,kx,ky,kz,ksquare);
dbx3=dbx3-nu*ksquare.*bx2;
dby3=dby3-nu*ksquare.*by2;
dbz3=dbz3-nu*ksquare.*bz2;
dbx3(2,1,1)=dbx3(2,1,1)-g/2;
dbx3(M,1,1)=dbx3(M,1,1)-g/2;
bx3=bx+deltat*dbx3;
by3=by+deltat*dby3;
bz3=bz+deltat*dbz3;

[dbx4,dby4,dbz4]=nonlinear(bx3,by3,bz3,kx,ky,kz,ksquare);
dbx4=dbx4-nu*ksquare.*bx3;
dby4=dby4-nu*ksquare.*by3;
dbz4=dbz4-nu*ksquare.*bz3;
dbx4(2,1,1)=dbx4(2,1,1)-g/2;
dbx4(M,1,1)=dbx4(M,1,1)-g/2;

bx=bx+deltat*(dbx1+2*dbx2+2*dbx3+dbx4)/6;
by=by+deltat*(dby1+2*dby2+2*dby3+dby4)/6;
bz=bz+deltat*(dbz1+2*dbz2+2*dbz3+dbz4)/6;


t=t+deltat


%  enforce exact conditions:
%  (1) vorticity is real
%	(2) its divergence is zero
bx=fft3(real(ifft3(bx)));
by=fft3(real(ifft3(by)));
bz=fft3(real(ifft3(bz)));
divg=kx.*bx+ky.*by+kz.*bz;
ksquare(1,1,1)=1;
bx=bx-kx.*divg./ksquare;
by=by-ky.*divg./ksquare;
bz=bz-kz.*divg./ksquare;
ksquare(1,1,1)=0;
%  end of enforcing code
end


end