I use gfortran, CodeBlocks environment
program main integer, parameter :: n = 3, nsteps = 100 real, parameter :: a = 0.0, b =10e-14 complex(8) :: x(0:n),h x = (/1.0, 0.0, 0.0,0.0/) h = (b - a)/nsteps print *,x call rk4sys(n,h,x,nsteps) end program main subroutine xpsys(n,x,f) complex(8), dimension (0:n) :: x complex(8), dimension (0:n) :: f integer n complex(8) :: E,d, hbar,W,omega0,omega,gamma1, delta d=3.335641E-30 hbar=(1.0546E-34,0.0) omega = 2.482/(6.24E18*hbar) W = 1.25E+10 omega0= 2.4743/(6.24E18*hbar) delta=0.0 gamma1= 1E15 E=(1/hbar)*W; f(0)= (0.,-1.)*d*(x(1)*E - x(2)*conjg(E))+gamma1*x(3) f(1)= conjg((0.,1.)*(x(2)*delta+ d*E*(x(0)- x(3)))-gamma1*x(2)/2); f(2)= (0., 1.)*(x(2)*delta + d*E*(x(0)- x(3)))-gamma1*x(2)/2 f(3)= (0., 1.)*d*(x(1)*E- x(2)*conjg(E))-gamma1*x(3) !DrhoDt(1)=-1i*d.*(rho(2).*(E) - rho(3).*conj(E))+gamma.*rho(4);+ !DrhoDt(2)=conj(1i.*(rho(3).*(delta)+ d.*E.*(rho(1)- rho(4)))-gamma.*rho(3)/2); !DrhoDt(3)=1i.*(rho(3).*(delta)+ d.*E.*(rho(1)- rho(4)))-gamma.*rho(3)/2; !DrhoDt(4)=1i.*d.*(rho(2).*(E)- rho(3).*conj(E))-gamma.*rho(4); end subroutine xpsys subroutine rk4sys(n,h,x,nsteps) complex(8), dimension (0:nsteps,0:3) :: Xmass complex(8) :: x(0:n) complex(8), allocatable :: y(:) complex(8), allocatable :: f(:,:) integer :: i, k, n complex(8) :: h allocate (y(0:n), f(0:n,4)) Do i=0,nsteps Do j=0,3 Xmass(i,j)=0 End Do End Do Xmass(0,0)=x(0) Xmass(0,1)=x(1) Xmass(0,2)=x(2) Xmass(0,3)=x(3) out: do k = 1,nsteps call xpsys(n,x,f(0,1)) in1: do i = 0,n y(i) = x(i) + 0.5*h*f(i,1) end do in1 call xpsys(n,y,f(0,2)) in2: do i = 0,n y(i) = x(i) + 0.5*h*f(i,2) end do in2 call xpsys(n,y,f(0,3)) in3: do i = 0,n y(i) = x(i) + h*f(i,3) end do in3 call xpsys(n,y,f(0,4)) in4: do i = 0,n x(i) = x(i) + (h/6.0)* (f(i,1) + 2.0*(f(i,2) + f(i,3)) + f(i,4)) Xmass(k,i) = x(i) end do in4 end do out k=0 OPEN(10,FILE='res.txt') !WRITE(10,'(5A25)') 'T','X1','X2','X3',' X4' DO k=0,nsteps print *,k, real (Xmass(k,0)),AIMAG(Xmass(k,0)), real (Xmass(k,1)), & AIMAG(Xmass(k,1)),real (Xmass(k,2)),AIMAG(Xmass(k,2)), & real (Xmass(k,3)),AIMAG(Xmass(k,3)) WRITE(10,'(I3, 8Z19.18)')k, real (Xmass(k,0)),AIMAG(Xmass(k,0)), real (Xmass(k,1)), & AIMAG(Xmass(k,1)),real (Xmass(k,2)),AIMAG(Xmass(k,2)), & real (Xmass(k,3)),AIMAG(Xmass(k,3)) END DO CLOSE(10) Xmass=0 end subroutine rk4sys This is displayed in the terminal - everything is in order.
1 0.90719088395362735 0.0000000000000000 0.00000000000000 -0.29033970820889698 0.0000000000000000 0.29033970820889698 9.2809116046372611E-002 0.0000000000000000000000000000001
In .txt is displayed - not the order
1 003FED07B52D39ECD6 00000000000000000000 000000000000000000 00BFD294ECFFDF317C 000000000000000000 003FD294ECFFDF317C 003FB7C2569630994D 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000BBDD294ECFFDF317C 000000000000000000003 ...