subroutine RKSYST(t0)	!, x0, xEnd)
use work
! 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)
real(8) t0
!real(8) x0(n), xEnd(n)

! Primele estimate ale lui DELTA-X:
!call DERIVS(x0,t0)
call DERIVS(t0)
w1 =h*f
xEnd =x0 +w1/2.D0
!do l=1,n
!xEnd(l) =x0(l) +w1(l)/2.
!enddo

! A 2-a estimare - in vectorul X:
t05 =t0+h/2._8
!call DERIVS(xEnd, t05)
call DERIVS(t05)
w2 =h*f
xEnd =x0 +w2/2.D0
!do l=1,n
!xEnd(l) =x0(l) +w2(l)/2.
!enddo

! A 3-a estimare:
!call DERIVS(xEnd, t05)
call DERIVS(t05)
w3 =h*f
xEnd =x0 +w3
!do l=1,n
!xEnd(l) =x0(l) +w3(l)
!enddo

! Ultima estimare:
!call DERIVS(xEnd, t0+h)
call DERIVS(t0+h)
w4 =h*f


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

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

end subroutine RKSYST
