#Arkaitz Rodas Bilbao
#© Derechos reservados excepto normas de uso y distribución de la UPV-EHU.

#---------------------------------------------------------------------------------------
#En este modulo realizamos la integracion numerica que nos permite obtener los 
#terminos del metodo runge-kutta
#---------------------------------------------------------------------------------------

disp('');

function cof=fc(HT,H1,d,db,tl,t,to,n,nt,hl)

 tu=to;
 
 hl=(t-to)/100;
 
 for p=1: 100
  
#--------------------------------------------------
#Interpolamos el hamiltoniano con una subrutina.
#--------------------------------------------------

 H2=fi(HT,tl,tu,n,nt);
 
#-----------------------------------------------------------
#Calculamos el producto tensorial de las matrices densidad.
#-----------------------------------------------------------
 
  D=f(d,db);
  
#----------------------------------------
#Obtenemos los terminos del integrando.
#----------------------------------------
  
  R=H2*D-D*H2;
  RT=H1*R-R*H1;

#-----------------------------
#Realizamos la traza parcial.
#-----------------------------
	
 Ru=fu(RT,n);
	
#---------------------------------------
#Finalmente integramos los terminos.
#---------------------------------------
	
  y1(p)=Ru(1,1);
  y2(p)=Ru(1,2);
  y3(p)=Ru(2,1);
  y4(p)=Ru(2,2);
  
   tu=tu+hl;
  
 endfor;
 
 k(1,1)=trapz(y1);
 k(1,2)=trapz(y2);
 k(2,1)=trapz(y3);
 k(2,2)=trapz(y4);
 
 cof=k;

endfunction

disp('');