program rungekutta
implicit none
integer, parameter :: dp = selected_real_kind(15,300)
integer :: i
real(kind=dp) z, y, t
do i=1,100
call rk2(dydt, dzdt, z, y, t)
print *, z, y
end do

contains
subroutine rk2(z, y, t)
implicit none
integer, parameter :: dp = selected_real_kind(15,300)
real(kind=dp) :: k1y, k1z, k2y, k2z, n, h=0.1, y=1, z=0, t
k1y = dydt(y,z,t)*h
k1z = dzdt(y,z,t)*h
k2z = dzdt(y + (0.5*k1y), z + (0.5*k1z), t + (0.5*h))*h
k2y = dydt(y, z +(0.5*k1z), t)*h
y = y + k2y
z = z + k2z
end subroutine

function dzdt(y,z,t)
real(kind=dp) :: y, z, t, dzdt, omega=1, A=0, B, C=0, D=0
B = omega**2
dzdt = A*y**3 + B*y - C*z + D*sin(omega*t)
end function

function dydt(z)
real(kind=dp) ::  z, dydt
dydt = z
end function
end program