module SigmoidPDE
   implicit none
   private
   public :: solve_sigmoid_pde

contains

   ! solve_sigmoid_pde
   !   Approximate solution of a sigmoid-based PDE in a hyperdimensional space
   !
   ! Parameters:
   !   u_in    - real array of dimension n: input vector (state at t=0)
   !   u_out   - real array of dimension n: output vector (state at T)
   !   dims    - integer: dimensionality of the space
   !   steps   - integer: number of time integration steps
   !   dt      - real: time step size
   !   epsilon - real: diffusion coefficient
   !
   ! Returns:
   !   u_out   - approximated solution after integration
   subroutine solve_sigmoid_pde(u_in, u_out, dims, steps, dt, epsilon)
      real, intent(in)    :: u_in(:)
      real, intent(out)   :: u_out(:)
      integer, intent(in) :: dims, steps
      real, intent(in)    :: dt, epsilon

      real, allocatable :: u(:), u_new(:)
      integer :: i, t

      ! Check input vector size
      if (size(u_in) /= dims) then
         print *, "Error: input vector size mismatch"
         stop
      end if

      allocate(u(dims), u_new(dims))
      u = u_in

      do t = 1, steps
         ! Compute next state via explicit Euler with diffusion + sigmoid reaction
         do i = 1, dims
            u_new(i) = u(i) + dt * (epsilon * laplacian(u, i, dims) + dsigmoid(u(i)))
         end do
         u = u_new
      end do

      u_out = u
      deallocate(u, u_new)

      return
   end subroutine solve_sigmoid_pde

   ! sigmoid
   !   Standard logistic function
   pure function sigmoid(x) result(s)
      real, intent(in) :: x
      real :: s

      s = 1.0 / (1.0 + exp(-x))
   end function sigmoid

   ! dsigmoid
   !   Derivative of the logistic function
   pure function dsigmoid(x) result(ds)
      real, intent(in) :: x
      real :: s
      real :: ds

      s  = sigmoid(x)
      ds = s * (1.0 - s)
   end function dsigmoid

   ! laplacian
   !   Approximates second spatial derivative for 1D index within hyperdimensional grid
   !   Here simplified as finite-difference with periodic boundary
   pure function laplacian(u, idx, n) result(lap)
      real, intent(in) :: u(:)
      integer, intent(in) :: idx, n
      integer :: im1, ip1
      real :: lap

      im1 = idx - 1
      ip1 = idx + 1
      if (im1 < 1) im1 = n
      if (ip1 > n) ip1 = 1

      lap = u(im1) - 2.0 * u(idx) + u(ip1)
   end function laplacian

end module SigmoidPDE

program main
   use SigmoidPDE
   implicit none

   integer, parameter :: dims = 100
   integer, parameter :: steps = 1000
   real,    parameter :: dt = 0.01
   real,    parameter :: epsilon = 0.1

   real :: u0(dims), uT(dims)
   integer :: i

   ! Initialize input vector
   do i = 1, dims
      u0(i) = 0.1 * i
   end do

   ! Solve PDE
   call solve_sigmoid_pde(u0, uT, dims, steps, dt, epsilon)

   ! Output first 10 elements of solution
   print *, 'Solution (first 10):'
   do i = 1, 10
      print '(F8.4)', uT(i)
   end do

end program main