function [tmat,SLin]=hbar_mrk4_quantum_linear_entropy(x1,p1,x2,E,lambda,zeta1,zeta2,N,hbar,tmax,dt,skipstep)
%close all
%clear all
format long %long %
%N=30
%lambda=0.0075
a=spdiags(sqrt(0:N-1)',1,N,N);
%%%%%%%%%%%%%%%%%%%%%Defining the Operators%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%id1=speye(N);
H0=hbar*(a'*a+0.5*speye(N,N));
H12=(hbar/2)*sqrt(lambda)*(a'+a)*(a'+a);
%%%%%%%%%%%%%%%%creating the state vector%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%psi_t=zeros(N,N);
%%%Check the function to learn to make a tensor product state vector%%%%%%
psi_t=hbar_sq_state_generator(x1,p1,x2,E,lambda,zeta1,zeta2,N,hbar);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Applying the Hamiltonian on the state vector%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%This tecnique is possible only for real symmetric Hamiltonian matrices%
%**************************************************************************
%**************************************************************************
%dt=0.0001
%tmax=2
skipstep=ceil(skipstep/(hbar));
tlen=tmax/(hbar*dt);
JJ=0
for indxt=0:tlen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (mod(skipstep+indxt,skipstep)==0) %&&(indxt*dt>=1.0)
%if (mod(10+indxt,10)==0) %just to test the program
rhoA=psi_t*psi_t';
dden=eig(full(rhoA));
%vind=find(abs(dden)>=10^-6);
%log_en=zeros(size(dden));
%log_en(vind)=log(abs(dden(vind)));
%rem_num=tlen-indxt
JJ=JJ+1
%svn(JJ)=-real(dden'*log_en)/log(2);
SLin(JJ)=1-sum(abs(dden).^2);
tmat(JJ)=indxt*dt*hbar;
%PR(JJ)=1/sum(abs(psi_t(:)).^4);
vn(JJ)=norm(psi_t(:));
%Hpsi=H0*psi_t+psi_t*H0+H12*(psi_t*H12);
%QEner(JJ)=psi_t(:)'*Hpsi(:)
%norm=sum(sum(abs(psi_t.^2)))
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if (mod(indxt,4*skipstep)==0)
% pointime=tmat(JJ)
% pointsvn=svn(JJ)
% Pratio=PR(JJ)
% Quantum_Energy=QEner(JJ)
% vector_norm=vn(JJ)
% rem_num=tlen-indxt
% progress=indxt/tlen*100
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%RK4th Method starts%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f1=sparse(-1i*dt*(H0*psi_t+psi_t*H0+H12*(psi_t*H12)));
f2=sparse(-1i*dt*(H0*(psi_t+0.5*f1)+(psi_t+0.5*f1)*H0+H12*((psi_t+0.5*f1)*H12)));
f3=sparse(-1i*dt*(H0*(psi_t+0.5*f2)+(psi_t+0.5*f2)*H0+H12*((psi_t+0.5*f2)*H12)));
f4=sparse(-1i*dt*(H0*(psi_t+f3)+(psi_t+f3)*H0+H12*((psi_t+f3)*H12)));
psi_t=sparse(psi_t+(1/6)*(f1+2*f2+2*f3+f4));
%psi_t=U*psi_t;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
%figure,plot(tmat,svn)
%save([num2str(fnum),'_',num2str(hbar),'_CupHam_Hbar_Ent'])
end