5 ! Contributed by Daniel Carrera
9 pure
subroutine rk4_vec(t
, Y
, dY
, h
)
10 real, intent(inout
) :: t
, Y(:)
12 real, dimension(size(Y
)) :: k1
, k2
, k3
, k4
15 pure
function dY(t0
, y0
)
17 real, intent(in
) :: t0
, y0(size(Y
))
23 k2
= dY(t
+ h
/2, Y
+ k1
*h
/2)
24 k3
= dY(t
+ h
/2, Y
+ k2
*h
/2)
25 k4
= dY(t
+ h
, Y
+ k3
*h
)
27 Y
= Y
+ (k1
+ 2*k2
+ 2*k3
+ k4
) * h
/6