docs: def write_rk4(path,name,laser,omega,gamma,r,Lij,states=None,verbose=1):
  	r"""
      This function writes the Fortran code needed to calculate the time evolution of the density matrix elements
      `\rho_{ij}` using the Runge-Kutta method of order 4.
  
      INPUT:
  
      - ``path`` - A string with the working directory where all files will be stored. It must end with ``/``.
  
      - ``name`` - A string with the name of the experiment. All files produced will begin with this name.
      
      - ``laser`` - A list of Laser objects (see the Laser class).
      
      - ``omega`` - A matrix or list of lists containing the frequency differences `\omega_{ij}`.
      
      - ``gamma`` - A matrix or list of lists containing the spontaneous decay frequencies `\gamma_{ij}`.
      
      - ``r`` - A list of three matrices or lists of lists containing the components of the position operator `r_{-1ij},r_{0ij},r_{1ij}`.
      
      - ``Lij`` - A list with elements of the form ``[i,j,[l1,l2,...]]`` representing the sets `L_{ij}` of which lasers excite wich transitions. It does not need to contain an element for all ``i,j`` pairs, but only those which have a laser that excites them.
      
      - ``Omega`` - A floating point number indicating the frequency scale for the equations. The frequencies ``omega`` and ``gamma`` are divided by this number. If ``None`` the equations and the input are taken in SI units.
  
      OUTPUT:
      
      - A file ``name.f90`` is created in ``path``.
      """
  
  	global omega_rescaled
  	
  	t0=time()
  	Ne=len(omega[0])
  	Nl=len(laser)
  
  	if states==None: states=range(1,Ne+1)
  
  	#We make some checks
  	for i in range(Ne):
  		for j in range(Ne):
  			b1=not ('.' in str(omega[i][j]) or 'e' in str(omega[i][j]))
  			if b1: raise ValueError,'omega must be composed of floating point numbers.'
  			b2=not ('.' in str(gamma[i][j]) or 'e' in str(gamma[i][j]))
  			if b2: raise ValueError,'gamma must be composed of floating point numbers.'
  
  	#We rescale the frequencies as requested.
  	#~ if Omega != None:
  		#~ omega_rescaled=[[omega[i][j]/Omega for j in range(Ne)] for i in range(Ne)]
  		#~ #gamma=[[gamma[i][j]/Omega for j in range(Ne)] for i in range(Ne)]
  	#~ else:
  		#~ omega_rescaled=omega[:]
  	omega_rescaled=omega[:]
  		
  
  	#We determine wether it is possible to eliminate explicit time-dependance
  	theta=find_phase_transformation(Ne,Nl,r,Lij)
  	#We find the detunings if required
  
  	#We construct the correspondence i <-> I between degenerate and non-degenerate indices.
  	i_d,I_nd,Nnd=calculate_iI_correspondence(omega)
  	#We get wich transitions each laser induces
  	detunings,detuningsij=laser_detunings(Lij,Nl,i_d,I_nd,Nnd)
  	
  	#We get how many transitions each  laser induces
  	detuning_indices=[len(detunings[i]) for i in range(Nl)]
  	
  	#The number of detunings
  	Nd=sum([len(detunings[l]) for l in range(Nl)])
  	combinations=detuning_combinations(detuning_indices)	
  
  	code0='''program evolution_rk4
  	implicit none
  	complex*16, dimension('''+str(Ne*(Ne+1)/2-1)+') :: x,k1,k2,k3,k4\n'
  	code0+='''    real*8 :: dt,t,ddelta,delta,delta0
  	integer :: i,j,n,ldelta,ndelta,detuning_index,n_aprox,n_mod
  
  	logical :: print_steps,run_spectrum\n'''
  	code0+='    real*8, dimension('+str(Nl)+') :: E0,detuning_knob\n'
  	code0+='    real*8, dimension('+str(Nd)+') :: detuning\n\n'
  
  	code0+="    open(unit=1,file='"+path+name+".dat',status='unknown')\n\n"
  	code0+='    n_aprox=1500\n'
  	code0+='    !We load the parameters\n'
  	code0+="    open(unit=2,file='"+path+name+"_params.dat',status='unknown')\n" 
  	code0+='''    read(2,*) n
      read(2,*) dt
      read(2,*) print_steps
      read(2,*) x
      read(2,*) E0\n'''
  
  	code0+='    read(2,*) detuning_knob\n'
  	code0+='    read(2,*) run_spectrum\n\n'
  	code0+='''    if (run_spectrum) then
  		read(2,*) ldelta
  		read(2,*) ndelta
  		read(2,*) ddelta
  		close(2)
  		delta0=detuning_knob(ldelta)
  		n_mod=ndelta/n_aprox
      else
  		ldelta=1; ndelta=1; ddelta=0; delta=0
  		close(2)
  		n_mod=n/n_aprox
      end if
      if (n_mod==0) n_mod=1\n\n\n'''
  
  
  	#We add the code to caculate all the initial detunings for each laser.
  	code0+='	!We calculate the initial detunings.\n'
  	#We find the minimal frequency corresponding to each laser.
  	omega_min,omega_min_indices=find_omega_min(omega_rescaled,Nl,detuningsij,i_d,I_nd)
  	
  	det_index=1
  	for l in range(Nl):
  		omega0=omega_min[l]
  		for p in detuningsij[l]:
  			code0+='	detuning('+str(det_index)+')='
  			code0+=format_double(omega0-omega_rescaled[p[0]][p[1]])+'+detuning_knob('+str(l+1)+')\n'
  			det_index+=1
  	code0+='\n'
  
  	code0+='''	t=0
  	if (.not. run_spectrum) WRITE(1,*) t,real(x),imag(x('''+str(Ne)+''':))\n'''
  	code0+='''	!We start the detuning variation\n'''
  	code0+='	delta=detuning_knob(ldelta)\n'
  	code0+='''    do j=1,ndelta
  		!We run the Runge Kutta method
  		t=0.0
  		do i=1,n-1\n'''
  
  	code0+='            call f(x          , t       , k1,   E0, detuning, detuning_knob)\n'
  	code0+='            call f(x+0.5*k1*dt, t+dt*0.5, k2,   E0, detuning, detuning_knob)\n'
  	code0+='            call f(x+0.5*k2*dt, t+dt*0.5, k3,   E0, detuning, detuning_knob)\n'
  	code0+='            call f(x    +k3*dt, t+dt    , k4,   E0, detuning, detuning_knob)\n'
  
  	code0+='''			x= x+(k1+2*k2+2*k3+k4)*dt/6
  			if (print_steps.and. .not. run_spectrum) print*,'t=',t,'delta=',delta
  			t= t+ dt
  			
  			if (isnan(real(x(1)))) stop 1
  			if (.not. run_spectrum .and. mod(i,n_mod)==0) WRITE(1,*) t,real(x),imag(x('''+str(Ne)+''':))
  		end do
  		if (print_steps) print*, 'delta=',delta,'percentage=',100*(delta-delta0)/(ddelta*ndelta)
  		
  		!We recalculate the detunings
  		if (run_spectrum) then
  			if (mod(j,n_mod)==0) WRITE(1,*) delta,real(x),imag(x('''+str(Ne)+''':))
  			delta=delta+ddelta
  			detuning_knob(ldelta)=detuning_knob(ldelta)+ddelta\n'''
  	
  
  
  	#We add the code to caculate all detunings for each laser
  	#This way of assigining a global index ll to the detunings ammounts to
  	#   ll=   number_of_previous_detunings 
  	#       + number_of_detuning_ordered_by_row_and_from_left_to_right_column
  	#like this
  	#->
  	#-> ->
  	#-> -> ->
  	#for each l
  	
  	#We find the minimal frequency corresponding to each laser		
  	omega_min,omega_min_indices=find_omega_min(omega_rescaled,Nl,detuningsij,i_d,I_nd)
  	
  	det_index=1
  	for l in range(Nl):
  		omega0=omega_min[l]
  		code0+='			if (ldelta=='+str(l+1)+') then\n'
  		for p in detuningsij[l]:
  			code0+='				detuning('+str(det_index)+')=detuning('+str(det_index)+')'
  			#code0+='+('+str(omega0-omega_rescaled[p[0]][p[1]])+'+ddelta\n'
  			code0+='+ddelta\n'
  			det_index+=1
  		code0+='			end if\n'
  	
  	code0+='''		end if
  	
  	
  	end do
  	
      close(1)
      
  end program\n\n'''
  
  
  	code0+='subroutine f(x,t,y,    E0, detuning,detuning_knob)\n'
  	
  	code0+='''    implicit none
      real*8, intent(in) :: t\n'''
  	code0+='    complex*16, dimension('+str(Ne*(Ne+1)/2-1)+'), intent(in)  :: x\n'
  	code0+='    complex*16, dimension('+str(Ne*(Ne+1)/2-1)+'), intent(out) :: y\n'
  	code0+='    real*8, dimension('+str(Nl)+'), intent(in) :: E0,detuning_knob\n'
  	code0+='    real*8, dimension('+str(Nd)+'), intent(in) :: detuning\n\n'
  
  	code0+='    complex*16 :: I,fact,aux\n'
  	code0+='    real*8 :: rho11\n\n'
  	code0+='    I=(0,1D0)\n'
  	
  	#We establish the scaling of the equations
  	#~ if Omega==None:
  		#~ h  =1.054571726e-34; e=1.602176565e-19
  		#~ code0+='    fact=I*'+str(e/h)+'\n'
  	#~ else:
  		#~ #code0+='    fact=I*'+str(float(Omega/sqrt(2)))+'\n'
  		#~ code0+='    fact=I*'+str(float(1/sqrt(2)))+'\n'
  		#~ #code0+='    fact=I*'+str(float(1/(sqrt(2)*Omega)))+'\n'
  	code0+='    fact=I*'+format_double(float(1/sqrt(2)))+'\n'
  	
  	#We give the code to calculate rho11
  	code0+='    rho11=1\n'
  	for i in range(1,Ne):
  		code0+='    rho11=rho11 -x('+str(i)+')\n'
  	code0+='\n\n'
  	
  	
  	####################################################################
  	#We produce the code for the first order equations.
  	####################################################################
  	if len(theta)>0:
  		code=''
  		for mu in range(1,Ne*(Ne+1)/2):
  			i,j,s=IJ(mu,Ne)
  			#print 'ecuacion mu=',mu,',i,j=',i,j
  			eqmu='    y('+str(mu)+')= 0\n'
  			
  			####################################################################
  			#We add the terms associated with the effective hamiltonian
  			#other than those associated with the phase transformation.
  			
  			for k in range(1,Ne+1):
  				#Case 1
  				if k>=j:
  					for l in Lij[i-1][k-1]:
  						if k>i:
  							#print 'E0^',l,-1,'r',i,k,'rho',k,j,'case 1.1'
  							eqmu+=add_line(Ne,mu,'+',laser,l,-1,  r,i,k,  k,j)
  						elif k<i:
  							#print 'E0^',l, 1,'r',i,k,'rho',k,j,'case 1.2'
  							eqmu+=add_line(Ne,mu,'+',laser,l, 1,  r,i,k,  k,j)
  							
  				#Case 2
  				elif k<j:
  					for l in Lij[i-1][k-1]:
  						if k>i:
  							#print 'E0^',l,-1,'r',i,k,'rhoa',j,k,'case 2.1'
  							eqmu+=add_line(Ne,mu,'+',laser,l,-1,  r,i,k,  j,k,True)
  						elif k<i:
  							#print 'E0^',l, 1,'r',i,k,'rhoa',j,k,'case 2.2'
  							eqmu+=add_line(Ne,mu,'+',laser,l, 1,  r,i,k,  j,k,True)
  				#Case 3
  				if k<=i:
  					for l in Lij[k-1][j-1]:
  						if k<j:
  							#print 'E0^',l,-1,'r',k,j,'rho',i,k,'case 3.1'
  							eqmu+=add_line(Ne,mu,'-',laser,l,-1,  r,k,j,  i,k)
  						elif k>j:
  							#print 'E0^',l, 1,'r',k,j,'rho',i,k,'case 3.2'
  							eqmu+=add_line(Ne,mu,'-',laser,l, 1,  r,k,j,  i,k)
  				#Case 4
  				elif k>i:
  					for l in Lij[k-1][j-1]:
  						if k<j:
  							#print 'E0^',l,-1,'r',k,j,'rhoa',k,i,'case 4.1'
  							eqmu+=add_line(Ne,mu,'-',laser,l,-1,  r,k,j,  k,i,True)
  						elif k>j:
  							#print 'E0^',l, 1,'r',k,j,'rhoa',k,i,'case 4.2'
  							eqmu+=add_line(Ne,mu,'-',laser,l, 1,  r,k,j,  k,i,True)
  
  			eqmu+='    y('+str(mu)+')=y('+str(mu)+')*fact\n'
  			
  			####################################################################
  			#We add the terms associated with the phase transformation.
  
  			extra=Theta(i,j,theta,omega_rescaled,omega_min,detunings,detuningsij,
  						combinations,detuning_indices,Lij,i_d,I_nd,Nnd,
  						verbose=verbose,states=states)
  
  			if extra!='':
  				eqmu+='    y('+str(mu)+')=y('+str(mu)+') + I*('+extra+')*x('+str(mu)+')\n'
  
  			####################################################################
  			
  			#~ if i==j:
  				#~ for k in range(1,Ne+1):
  					#~ if k < i:
  						#~ muii=Mu(i,i,s=1,N=Ne)
  						#~ eqmu+='    y('+str(mu)+')= y('+str(mu)+') - ('+format_double(gamma[i-1][k-1])+')*x('+str(muii)+')\n'
  					#~ elif k > i:
  						#~ mukk=Mu(k,k,s=1,N=Ne)
  						#~ eqmu+='    y('+str(mu)+')= y('+str(mu)+') - ('+format_double(gamma[i-1][k-1])+')*x('+str(mukk)+')\n'
  				#~ eqmu+='\n'
  			#~ else:
  				#~ eqmu+='    y('+str(mu)+')= y('+str(mu)+') - ('+format_double(gamma[i-1][j-1]/2)+')*x('+str(mu)+')\n'
  #~ 
  
  			####################################################################
  			code+=eqmu+'\n'
  
  		#We add the terms associated with spontaneous decay.		
  		#First for populations.
  		for i in range(2,Ne+1):
  			mu=Mu(i,i,1,Ne)
  			for k in range(1,Ne+1):
  				gams=0
  				if k<i:
  					gams+=gamma[i-1][k-1]
  				elif k>i:
  					nu=Mu(k,k,1,Ne)
  					ga=gamma[i-1][k-1]
  					if ga != 0:
  						code+='    y('+str(mu)+')=y('+str(mu)+')'
  						code+='-('+format_double(ga)+')*x('+str(nu)+')\n'
  				if gams!=0:
  					code+='    y('+str(mu)+')=y('+str(mu)+')'
  					code+='-('+format_double(gams)+')*x('+str(mu)+')\n'
  
  		#And now for coherences	
  		for i in range(1,Ne+1):
  			for j in range(1,i):
  				gams=gamma[i-1][j-1]/2
  				if gams!=0:
  					for a in range(i+1,Ne+1):
  						mu=Mu(a,i,+1,Ne)
  						code+='    y('+str(mu)+')=y('+str(mu)+')'
  						code+='-('+format_double(gams)+')*x('+str(mu)+')\n'
  
  						#~ mu=Mu(a,i,-1,Ne)
  						#~ code+='    y('+str(mu)+')=y('+str(mu)+')'
  						#~ code+='-('+format_double(gams)+')*x('+str(mu)+')\n'
  
  					for b in range(1,i):
  						mu=Mu(i,b,+1,Ne)
  						code+='    y('+str(mu)+')=y('+str(mu)+')'
  						code+='-('+format_double(gams)+')*x('+str(mu)+')\n'
  
  						#~ mu=Mu(i,b,-1,Ne)
  						#~ code+='    y('+str(mu)+')=y('+str(mu)+')'
  						#~ code+='-('+format_double(gams)+')*x('+str(mu)+')\n'
  
  
  
  
  		
  		####################################################################
  		####################################################################
  		####################################################################
  		#code+='    y=y/'+str(Omega)+'\n'
  		
  		f=file(path+name+'.f90','w')
  		code=code0+code+'end subroutine\n'
  		f.write(code)
  		f.close()
  				
  		return time()-t0
  	else:
  		print 'There was no phase transformation capable of eliminating explicit time dependance.'