program pendulum_forward_euler implicit none double precision, parameter :: pi = acos(-1.0d0) double precision :: q(2) double precision :: h, end_time integer :: n, n_steps end_time = 10 * pi ! End time h = 0.1d0 ! Step size n_steps = nint(end_time / h) ! Number of steps q(:) = [0.01d0 * pi, 0.0d0] ! Initial state q = (theta, phi) ! Forward Euler do n = 1, n_steps q = q + h * [q(2), -sin(q(1))] ! print time, theta print *, n * h, q(1) end do end program pendulum_forward_euler