subroutine RKSYST(t0)
use work
! DOUBLE
! Runge-Kutta pentru sistem de ordinul "n"
! dx1/dt =f1(x,t); dx2/dt =f2(x,t); ...
! * DERIVS = subrutina care calculeaza derivatele; se declara EXTERNAL.
! * call DERIVS(x,t,f,n)
! t0: t-initial
! h: pasul
! x: tablou pentru valorile initiale x(n)
! xEnd:tablou pentru valorile finale x(n)
! w: tablou pentru valorile intermediare w(4,n)
! n: numarul ecuatiilor
! F: tablou pentru valorile derivatelor

IMPLICIT REAL*8 (a-h,o-z)

! First estimates of DELTA-X:
call DERIVS(t0)
w1 =h*f
!xEnd =x0 +w1/2.
do l=1,n
xEnd(l) =x0(l) +w1(l)/2.
enddo

! 2nd estimate - in vector X:
t05 =t0+h/2.
call DERIVS(t05)	! +h/2.)
w2 =h*f
!xEnd =x0 +w2/2.
do l=1,n
xEnd(l) =x0(l) +w2(l)/2.
enddo

! 3rd estimate:
call DERIVS(t05)	! +h/2.)
w3 =h*f
!xEnd =x0 +w3
do l=1,n
xEnd(l) =x0(l) +w3(l)
enddo

! Last estimate:
call DERIVS(t0+h)
w4 =h*f

do i=1,n
xEnd(i) =x0(i) +( w1(i) +2.*(w2(i) +w3(i)) +w4(i) )/6.
enddo

!xEnd =x0 +( w1 +2.*(w2 +w3) +w4 )/6.

end subroutine RKSYST
