program logistic_sigmoid_pde
	implicit none
	integer, parameter :: NDIM = 4       ! 次元数
	integer, parameter :: N = 50         ! 各次元の格子数
	integer :: d, idx, idxp, idxm
	integer, dimension(NDIM) :: sizes, strides
	integer :: total_points
	real(kind=8), allocatable :: u(:), u_new(:)
	real(kind=8) :: D, r, dt, dx, t, t_end

	! 格子情報の初期化
	total_points = 1
	do d = 1, NDIM
		sizes(d) = N
		total_points = total_points * N
	end do

	allocate(u(total_points), u_new(total_points))

	! フラット配列用ストライド計算
	strides(1) = 1
	do d = 2, NDIM
		strides(d) = strides(d-1) * sizes(d-1)
	end do

	! 物理パラメータ設定
	dx    = 1.0d0 / (N - 1)
	D     = 0.1d0                       ! 拡散係数
	r     = 1.0d0                       ! 反応（ロジスティック）係数
	dt    = 0.1d0 * dx * dx / D         ! 安定性条件に基づく時間ステップ
	t_end = 1.0d0                       ! シミュレーション終了時間

	! 初期条件: 乱数
	call random_seed()
	call random_number(u)

	t = 0.0d0
	do while (t < t_end)
		! 時間発展（陽解法）
		do idx = 1, total_points
			! 反応項：u*(1-u)
			u_new(idx) = r * u(idx) * (1.0d0 - u(idx))
			! 拡散項：ラプラシアンの有限差分近似
			do d = 1, NDIM
				idxp = idx + strides(d)
				idxm = idx - strides(d)
				! 周期境界条件
				if (mod(idx-1, strides(d)*sizes(d)) == 0) then
					idxm = idx + strides(d)*(sizes(d)-1)
				end if
				if (mod(idx-1, strides(d)*sizes(d)) >= strides(d)*(sizes(d)-1)) then
					idxp = idx - strides(d)*(sizes(d)-1)
				end if
				u_new(idx) = u_new(idx) + D * &
				              (u(idxp) - 2.0d0*u(idx) + u(idxm)) / (dx*dx)
			end do
			! 次ステップへの更新
			u_new(idx) = u(idx) + dt * u_new(idx)
		end do

		u = u_new
		t = t + dt
	end do

	! 結果出力（例：先頭10点）
	print *, 't = ', t
	do idx = 1, 10
		print *, idx, u(idx)
	end do

	deallocate(u, u_new)
end program logistic_sigmoid_pde