diff --git a/Analysis/maxent_wrapper_mod.F90 b/Analysis/maxent_wrapper_mod.F90
index 09e499fa..178452c9 100644
--- a/Analysis/maxent_wrapper_mod.F90
+++ b/Analysis/maxent_wrapper_mod.F90
@@ -43,7 +43,13 @@ Real (Kind=Kind(0.d0)) function XKER_ph(tau,om, beta)
 
        pi = 3.1415927
 
-       XKER_ph = (exp(-tau*om) + exp(-( beta - tau )*om ) )/( pi*(1.d0 + exp( - beta * om ) ) )
+       if ( beta*om > 600 ) then
+         XKER_ph = (exp(-tau*om) + exp(-( beta - tau )*om )) / pi
+       elseif ( beta*om < -600 ) then
+         XKER_ph = (exp((beta-tau)*om) + exp(tau*om)) / pi
+       else
+         XKER_ph = (exp(-tau*om) + exp(-( beta - tau )*om ) )/( pi*(1.d0 + exp( - beta * om ) ) )
+       endif
 
      end function XKER_ph
      
@@ -54,7 +60,13 @@ Real (Kind=Kind(0.d0)) function XKER_ph_c(tau,om, beta)
 
        pi = 3.1415927
 
-       XKER_ph_c = (exp(-tau*om) + exp(-( beta - tau )*om ) )/( pi*(1.d0 + exp( - beta * om ) ) )
+       if ( beta*om > 600 ) then
+         XKER_ph_c = (exp(-tau*om) + exp(-( beta - tau )*om )) / pi
+       elseif ( beta*om < -600 ) then
+         XKER_ph_c = (exp((beta-tau)*om) + exp(tau*om)) / pi
+       else
+         XKER_ph_c = (exp(-tau*om) + exp(-( beta - tau )*om ) )/( pi*(1.d0 + exp( - beta * om ) ) )
+       endif
 
      end function XKER_ph_c
 
@@ -65,7 +77,13 @@ Real (Kind=Kind(0.d0)) function XKER_pp(tau,om, beta)
 
        pi = 3.1415927
 
-       XKER_pp = exp(-tau*om) / ( pi*(1.d0 + exp( - beta * om ) ) )
+       if ( beta*om > 600 ) then
+         XKER_pp = exp(-tau*om) / pi
+       elseif ( beta*om < -600 ) then
+         XKER_pp = exp((beta-tau)*om) / pi
+       else
+         XKER_pp = exp(-tau*om) / ( pi*(1.d0 + exp( - beta * om ) ) )
+       endif
 
      end function XKER_pp
 
@@ -76,7 +94,13 @@ Real (Kind=Kind(0.d0)) function XKER_p(tau,om, beta)
 
        pi = 3.1415927
 
-       XKER_p  = exp(-tau*om) / ( pi*(1.d0 + exp( - beta * om ) ) )
+       if ( beta*om > 600 ) then
+         XKER_p  = exp(-tau*om) / pi
+       elseif ( beta*om < -600 ) then
+         XKER_p  = exp( (beta-tau)*om) / pi
+       else
+         XKER_p  = exp(-tau*om) / ( pi*(1.d0 + exp( - beta * om ) ) )
+       endif
 
      end function XKER_p
 
@@ -102,7 +126,11 @@ Real (Kind=Kind(0.d0)) function F_QFI_ph(om, beta)
       real (Kind=Kind(0.d0)) ::  om, beta
       real (Kind=Kind(0.d0)) :: pi
       pi = 3.1415927
-      F_QFI_ph = (4.d0/pi) * ( (exp(beta*om) - 1.d0)/( exp(beta*om) + 1.d0 ) )**2
+      if ( abs(beta*om) > 600 ) then
+        F_QFI_ph = 4.d0/pi
+      else
+        F_QFI_ph = (4.d0/pi) * ( (exp(beta*om) - 1.d0)/( exp(beta*om) + 1.d0 ) )**2
+      endif
 
      end function F_QFI_ph
      
@@ -112,7 +140,11 @@ Real (Kind=Kind(0.d0)) function F_QFI_ph_c(om, beta)
       real (Kind=Kind(0.d0)) ::  om, beta
       real (Kind=Kind(0.d0)) :: pi
       pi = 3.1415927
-      F_QFI_ph_c = (4.d0/pi) * ( (exp(beta*om) - 1.d0)/( exp(beta*om) + 1.d0 ) )**2
+      if ( abs(beta*om) > 600 ) then
+        F_QFI_ph_c = 4.d0/pi
+      else
+        F_QFI_ph_c = (4.d0/pi) * ( (exp(beta*om) - 1.d0)/( exp(beta*om) + 1.d0 ) )**2
+      endif
 
      end function F_QFI_ph_c
 
@@ -121,7 +153,11 @@ Real (Kind=Kind(0.d0)) function F_QFI_pp(om, beta)
       real (Kind=Kind(0.d0)) ::  om, beta
       real (Kind=Kind(0.d0)) :: pi
       pi = 3.1415927
-      F_QFI_pp = (4.d0/pi) * ( (exp(beta*om) - 1.d0)/( exp(beta*om) + 1.d0 ) )**2
+      if ( abs(beta*om) > 600 ) then
+        F_QFI_pp = 4.d0/pi
+      else
+        F_QFI_pp = (4.d0/pi) * ( (exp(beta*om) - 1.d0)/( exp(beta*om) + 1.d0 ) )**2
+      endif
 
      end function F_QFI_pp
 
@@ -262,7 +298,13 @@ Real (Kind=Kind(0.d0)) function Back_trans_ph(Aom, om, beta )
        Implicit None
        real (Kind=Kind(0.d0)), intent(in) ::  Aom, beta, om
 
-       Back_trans_ph = Aom/(1.d0 + exp(-beta*om) )
+       if ( beta*om > 600 ) then
+         Back_trans_ph = Aom
+       elseif ( beta*om < -600 ) then
+         Back_trans_ph = Aom * exp(beta*om)
+       else
+         Back_trans_ph = Aom/(1.d0 + exp(-beta*om) )
+       endif
        ! This gives S(q,om) = chi(q,om)/(1 - e^(-beta om))
 
      end function BACK_TRANS_PH
@@ -276,6 +318,10 @@ Real (Kind=Kind(0.d0)) function Back_trans_ph_c(Aom, om, beta)
        Zero = 1.D-8
        if ( abs(om) < zero ) then
           Back_trans_ph_c = beta * Aom/2.d0
+       elseif ( beta*om > 600 ) then
+          Back_trans_ph_c = Aom / om
+       elseif ( beta*om < -600 ) then
+          Back_trans_ph_c = -Aom / om
        else
           Back_trans_ph_c = Aom * (1.d0 - exp(-beta*om) ) / (om *( 1.d0 + exp(-beta*om) ) )
           ! This gives sigma'(q,om) = A_c(q,om)*(1 - e^(-beta om))/(1 + e^(-beta om))/om
@@ -296,6 +342,10 @@ Real (Kind=Kind(0.d0)) function Back_trans_pp(Aom, om, beta)
        Zero = 1.D-8
        if ( abs(om) < Zero ) then
           Back_trans_pp = beta * Aom/2.d0
+       elseif ( beta*om > 600 ) then
+          Back_trans_pp = Aom / om
+       elseif ( beta*om < -600 ) then
+          Back_trans_pp = -Aom / om
        else
           Back_trans_pp = Aom * (1.d0 - exp(-beta*om) ) / (om *( 1.d0 + exp(-beta*om) ) )
        endif
@@ -313,7 +363,14 @@ Real (Kind=Kind(0.d0)) function XKER_p_ph(tau,om, beta)
        real (Kind=Kind(0.d0)) :: tau, om, pi, beta
 
        pi = 3.1415927
-       XKER_p_ph  =  (exp(-tau*om)  + exp(-(beta-tau)*om)) / (pi*(1.d0 + exp( -beta * om )) )
+
+       if ( beta*om > 600 ) then
+         XKER_p_ph = (exp(-tau*om) + exp(-(beta-tau)*om)) / pi
+       elseif ( beta*om < -600 ) then
+         XKER_p_ph = (exp((beta-tau)*om) + exp(tau*om)) / pi
+       else
+         XKER_p_ph = (exp(-tau*om) + exp(-(beta-tau)*om)) / (pi*(1.d0 + exp( -beta * om )) )
+       endif
 
      end function XKER_p_ph
 
@@ -442,10 +499,16 @@ Subroutine Set_default(Default,beta,Channel, OM_st, Om_en, xmom1,Default_model_e
        case("PH")
          If (.not. Default_model_exists ) Default = 1.d0/(Om_en - Om_st) ! Flat  default   
          !Compute   sum rule  for  A(om)
+         ! Note: PH channel uses symmetric spectrum with Om_st = 0 (enforced in Max_SAC.F90),
+         ! so om > 0 always holds here and the beta*om < -600 branch is never reached.
          X  = 0.d0
          Do  nw = 1, Ndis
              Om = Om_st + dble(nw)*dom
-             Default(nw)  = (1.d0 + exp(-beta*om)) * Default(nw)
+             if ( beta*om > 600 ) then
+               Default(nw) = Default(nw)
+             else
+               Default(nw) = (1.d0 + exp(-beta*om)) * Default(nw)
+             endif
              X = X + Default(nw) 
          enddo
          X = X*dom
@@ -467,6 +530,10 @@ Subroutine Set_default(Default,beta,Channel, OM_st, Om_en, xmom1,Default_model_e
                Om = Om_st + dble(nw)*dom
                if ( abs(om) < zero ) then
                   Default(nw) = Default(nw)*2.d0/ beta 
+               elseif ( beta*om > 600 ) then
+                  Default(nw) = Default(nw) * om
+               elseif ( beta*om < -600 ) then
+                  Default(nw) = Default(nw) * (-om)
                else
                   Default(nw) = Default(nw) * (om *( 1.d0 + exp(-beta*om) ) )/ (1.d0 - exp(-beta*om) ) 
                endif
@@ -492,6 +559,10 @@ Subroutine Set_default(Default,beta,Channel, OM_st, Om_en, xmom1,Default_model_e
                Om = Om_st + dble(nw)*dom
                if ( abs(om) < zero ) then
                   Default(nw) = Default(nw)*2.d0/ beta 
+               elseif ( beta*om > 600 ) then
+                  Default(nw) = Default(nw) * om
+               elseif ( beta*om < -600 ) then
+                  Default(nw) = Default(nw) * (-om)
                else
                   Default(nw) = Default(nw) * (om *( 1.d0 + exp(-beta*om) ) )/ (1.d0 - exp(-beta*om) ) 
                endif
diff --git a/testsuite/Analysis.tests/1-maxent-kernel-overflow.F90 b/testsuite/Analysis.tests/1-maxent-kernel-overflow.F90
new file mode 100644
index 00000000..8437e33a
--- /dev/null
+++ b/testsuite/Analysis.tests/1-maxent-kernel-overflow.F90
@@ -0,0 +1,292 @@
+! Test that the overflow-protected kernel and helper functions in
+! MaxEnt_Wrapper_mod are smooth across the beta*omega = +-600 branch
+! threshold and produce finite, correct results at extreme arguments.
+!
+! Strategy
+! --------
+!  (a) Continuity: at beta*omega = +-599 (just inside the exact branch)
+!      the exact formula and the asymptotic formula must agree to machine
+!      precision, because the neglected terms are O(exp(-599)) ~ 1e-260.
+!      We compute both and require relative agreement < 1e-12.
+!
+!  (b) Extreme values: at beta*omega = +-1000 the naive formula would
+!      overflow.  We verify the function returns a finite value matching
+!      the hand-coded asymptotic result.
+
+Program Test_MaxEnt_Kernel_Overflow
+
+  Use MaxEnt_Wrapper_mod
+  implicit none
+
+  real(Kind=Kind(0.d0)) :: beta, tau, om, pi
+  real(Kind=Kind(0.d0)) :: val, ref, relerr
+  real(Kind=Kind(0.d0)) :: tol
+  real(Kind=Kind(0.d0)) :: taus(4)
+  integer :: it, nerr
+
+  ! Use the same single-precision pi constant as maxent_wrapper_mod
+  pi  = 3.1415927
+  beta = 10.d0
+  tol  = 1.d-12
+  nerr = 0
+
+  taus(1) = 0.d0
+  taus(2) = beta/4.d0
+  taus(3) = beta/2.d0
+  taus(4) = beta
+
+  ! =====================================================================
+  ! XKER_ph  (and implicitly XKER_ph_c which is the same expression)
+  ! Exact: (exp(-tau*om) + exp(-(beta-tau)*om)) / (pi*(1+exp(-beta*om)))
+  ! beta*om >> 1 :  denom -> pi   =>  (exp(-tau*om)+exp(-(beta-tau)*om))/pi
+  ! beta*om << -1:  multiply by exp(beta*om) => (exp((beta-tau)*om)+exp(tau*om))/pi
+  ! =====================================================================
+  write(*,*) "--- XKER_ph ---"
+  do it = 1, 4
+     tau = taus(it)
+     ! (a) Continuity at +599
+     om  = 599.d0/beta
+     val = XKER_ph(tau, om, beta)
+     ref = (exp(-tau*om) + exp(-(beta-tau)*om)) / pi   ! asymptotic
+     call check("XKER_ph +599", val, ref, tol, nerr)
+
+     ! (a) Continuity at -599
+     om  = -599.d0/beta
+     val = XKER_ph(tau, om, beta)
+     ref = (exp((beta-tau)*om) + exp(tau*om)) / pi
+     call check("XKER_ph -599", val, ref, tol, nerr)
+
+     ! (b) Extreme +1000
+     om  = 1000.d0/beta
+     val = XKER_ph(tau, om, beta)
+     ref = (exp(-tau*om) + exp(-(beta-tau)*om)) / pi
+     call check("XKER_ph +1000", val, ref, tol, nerr)
+
+     ! (b) Extreme -1000
+     om  = -1000.d0/beta
+     val = XKER_ph(tau, om, beta)
+     ref = (exp((beta-tau)*om) + exp(tau*om)) / pi
+     call check("XKER_ph -1000", val, ref, tol, nerr)
+  enddo
+
+  ! =====================================================================
+  ! XKER_pp
+  ! Exact: exp(-tau*om) / (pi*(1+exp(-beta*om)))
+  ! beta*om >> 1 :  exp(-tau*om)/pi
+  ! beta*om << -1:  exp((beta-tau)*om)/pi
+  ! =====================================================================
+  write(*,*) "--- XKER_pp ---"
+  do it = 1, 4
+     tau = taus(it)
+     om  = 599.d0/beta
+     val = XKER_pp(tau, om, beta)
+     ref = exp(-tau*om) / pi
+     call check("XKER_pp +599", val, ref, tol, nerr)
+
+     om  = -599.d0/beta
+     val = XKER_pp(tau, om, beta)
+     ref = exp((beta-tau)*om) / pi
+     call check("XKER_pp -599", val, ref, tol, nerr)
+
+     om  = 1000.d0/beta
+     val = XKER_pp(tau, om, beta)
+     ref = exp(-tau*om) / pi
+     call check("XKER_pp +1000", val, ref, tol, nerr)
+
+     om  = -1000.d0/beta
+     val = XKER_pp(tau, om, beta)
+     ref = exp((beta-tau)*om) / pi
+     call check("XKER_pp -1000", val, ref, tol, nerr)
+  enddo
+
+  ! =====================================================================
+  ! XKER_p  (already had the fix; verify it still works)
+  ! =====================================================================
+  write(*,*) "--- XKER_p ---"
+  do it = 1, 4
+     tau = taus(it)
+     om  = 599.d0/beta
+     val = XKER_p(tau, om, beta)
+     ref = exp(-tau*om) / pi
+     call check("XKER_p +599", val, ref, tol, nerr)
+
+     om  = -599.d0/beta
+     val = XKER_p(tau, om, beta)
+     ref = exp((beta-tau)*om) / pi
+     call check("XKER_p -599", val, ref, tol, nerr)
+
+     om  = 1000.d0/beta
+     val = XKER_p(tau, om, beta)
+     ref = exp(-tau*om) / pi
+     call check("XKER_p +1000", val, ref, tol, nerr)
+
+     om  = -1000.d0/beta
+     val = XKER_p(tau, om, beta)
+     ref = exp((beta-tau)*om) / pi
+     call check("XKER_p -1000", val, ref, tol, nerr)
+  enddo
+
+  ! =====================================================================
+  ! XKER_p_ph  (same expression as XKER_ph)
+  ! =====================================================================
+  write(*,*) "--- XKER_p_ph ---"
+  do it = 1, 4
+     tau = taus(it)
+     om  = 599.d0/beta
+     val = XKER_p_ph(tau, om, beta)
+     ref = (exp(-tau*om) + exp(-(beta-tau)*om)) / pi
+     call check("XKER_p_ph +599", val, ref, tol, nerr)
+
+     om  = -599.d0/beta
+     val = XKER_p_ph(tau, om, beta)
+     ref = (exp((beta-tau)*om) + exp(tau*om)) / pi
+     call check("XKER_p_ph -599", val, ref, tol, nerr)
+
+     om  = 1000.d0/beta
+     val = XKER_p_ph(tau, om, beta)
+     ref = (exp(-tau*om) + exp(-(beta-tau)*om)) / pi
+     call check("XKER_p_ph +1000", val, ref, tol, nerr)
+
+     om  = -1000.d0/beta
+     val = XKER_p_ph(tau, om, beta)
+     ref = (exp((beta-tau)*om) + exp(tau*om)) / pi
+     call check("XKER_p_ph -1000", val, ref, tol, nerr)
+  enddo
+
+  ! =====================================================================
+  ! F_QFI_ph   (= (4/pi)*tanh(beta*om/2)^2 )
+  ! |beta*om| >> 1 :  tanh -> +-1,  squared -> 1  =>  4/pi
+  ! =====================================================================
+  write(*,*) "--- F_QFI_ph ---"
+  om  = 599.d0/beta
+  val = F_QFI_ph(om, beta)
+  ref = 4.d0/pi
+  call check("F_QFI_ph +599", val, ref, tol, nerr)
+
+  om  = -599.d0/beta
+  val = F_QFI_ph(om, beta)
+  ref = 4.d0/pi
+  call check("F_QFI_ph -599", val, ref, tol, nerr)
+
+  om  = 1000.d0/beta
+  val = F_QFI_ph(om, beta)
+  ref = 4.d0/pi
+  call check("F_QFI_ph +1000", val, ref, tol, nerr)
+
+  om  = -1000.d0/beta
+  val = F_QFI_ph(om, beta)
+  ref = 4.d0/pi
+  call check("F_QFI_ph -1000", val, ref, tol, nerr)
+
+  ! =====================================================================
+  ! Back_trans_ph   Aom/(1+exp(-beta*om))
+  ! beta*om >> 1 :  Aom
+  ! beta*om << -1:  Aom*exp(beta*om)
+  ! =====================================================================
+  write(*,*) "--- Back_trans_ph ---"
+  om  = 599.d0/beta
+  val = Back_trans_ph(1.d0, om, beta)
+  ref = 1.d0
+  call check("Back_trans_ph +599", val, ref, tol, nerr)
+
+  om  = -599.d0/beta
+  val = Back_trans_ph(1.d0, om, beta)
+  ref = exp(beta*om)
+  call check("Back_trans_ph -599", val, ref, tol, nerr)
+
+  om  = 1000.d0/beta
+  val = Back_trans_ph(1.d0, om, beta)
+  ref = 1.d0
+  call check("Back_trans_ph +1000", val, ref, tol, nerr)
+
+  om  = -1000.d0/beta
+  val = Back_trans_ph(1.d0, om, beta)
+  ref = exp(beta*om)
+  call check("Back_trans_ph -1000", val, ref, tol, nerr)
+
+  ! =====================================================================
+  ! Back_trans_ph_c   Aom*tanh(beta*om/2)/om
+  ! beta*om >> 1 :  Aom/om
+  ! beta*om << -1:  -Aom/om
+  ! =====================================================================
+  write(*,*) "--- Back_trans_ph_c ---"
+  om  = 599.d0/beta
+  val = Back_trans_ph_c(1.d0, om, beta)
+  ref = 1.d0/om
+  call check("Back_trans_ph_c +599", val, ref, tol, nerr)
+
+  om  = -599.d0/beta
+  val = Back_trans_ph_c(1.d0, om, beta)
+  ref = -1.d0/om
+  call check("Back_trans_ph_c -599", val, ref, tol, nerr)
+
+  om  = 1000.d0/beta
+  val = Back_trans_ph_c(1.d0, om, beta)
+  ref = 1.d0/om
+  call check("Back_trans_ph_c +1000", val, ref, tol, nerr)
+
+  om  = -1000.d0/beta
+  val = Back_trans_ph_c(1.d0, om, beta)
+  ref = -1.d0/om
+  call check("Back_trans_ph_c -1000", val, ref, tol, nerr)
+
+  ! =====================================================================
+  ! Back_trans_pp   same expression as Back_trans_ph_c
+  ! =====================================================================
+  write(*,*) "--- Back_trans_pp ---"
+  om  = 599.d0/beta
+  val = Back_trans_pp(1.d0, om, beta)
+  ref = 1.d0/om
+  call check("Back_trans_pp +599", val, ref, tol, nerr)
+
+  om  = -599.d0/beta
+  val = Back_trans_pp(1.d0, om, beta)
+  ref = -1.d0/om
+  call check("Back_trans_pp -599", val, ref, tol, nerr)
+
+  om  = 1000.d0/beta
+  val = Back_trans_pp(1.d0, om, beta)
+  ref = 1.d0/om
+  call check("Back_trans_pp +1000", val, ref, tol, nerr)
+
+  om  = -1000.d0/beta
+  val = Back_trans_pp(1.d0, om, beta)
+  ref = -1.d0/om
+  call check("Back_trans_pp -1000", val, ref, tol, nerr)
+
+  ! =====================================================================
+  ! Summary
+  ! =====================================================================
+  if (nerr > 0) then
+     write(*,'(A,I0,A)') "FAILED: ", nerr, " check(s) exceeded tolerance"
+     stop 2
+  endif
+  write(*,*) "success"
+
+contains
+
+  subroutine check(label, val, ref, tol, nerr)
+    use ieee_arithmetic, only: ieee_is_finite
+    character(len=*), intent(in) :: label
+    real(Kind=Kind(0.d0)), intent(in) :: val, ref, tol
+    integer, intent(inout) :: nerr
+    real(Kind=Kind(0.d0)) :: relerr, denom
+
+    ! Check finiteness before computing relerr to avoid NaN/Inf propagation
+    if (.not. ieee_is_finite(val)) then
+       write(*,'(A,A,A,ES14.6)') "  FAIL ", label, ": non-finite val=", val
+       nerr = nerr + 1
+       return
+    endif
+
+    denom = max(abs(ref), abs(val), tiny(1.d0))
+    relerr = abs(val - ref) / denom
+
+    if (relerr > tol) then
+       write(*,'(A,A,A,ES14.6,A,ES14.6,A,ES10.2)') &
+            "  FAIL ", label, ": val=", val, " ref=", ref, " relerr=", relerr
+       nerr = nerr + 1
+    endif
+  end subroutine check
+
+end Program Test_MaxEnt_Kernel_Overflow
diff --git a/testsuite/Analysis.tests/CMakeLists.txt b/testsuite/Analysis.tests/CMakeLists.txt
new file mode 100644
index 00000000..7060a0d7
--- /dev/null
+++ b/testsuite/Analysis.tests/CMakeLists.txt
@@ -0,0 +1,19 @@
+cmake_minimum_required(VERSION 3.5...4.2)
+Project(Analysis.tests C Fortran)
+
+add_executable(1-maxent-kernel-overflow 1-maxent-kernel-overflow.F90)
+
+set(analysistests 1-maxent-kernel-overflow)
+
+foreach(iter ${analysistests})
+  foreach(iflag ${FLAG_LIST})
+    target_compile_options(${iter} PRIVATE ${iflag})
+  endforeach(iflag)
+  target_link_libraries(${iter} ${CMAKE_CURRENT_SOURCE_DIR}/../../Analysis/maxent_wrapper_mod.o)
+  target_link_libraries(${iter} ${LIB_LIST})
+  target_include_directories(${iter} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../../Analysis/)
+  target_include_directories(${iter} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../../Libraries/Modules/)
+endforeach(iter)
+
+enable_testing()
+add_test(1-maxent-kernel-overflow 1-maxent-kernel-overflow)
diff --git a/testsuite/CMakeLists.txt b/testsuite/CMakeLists.txt
index 6a8ea471..521bc4ad 100644
--- a/testsuite/CMakeLists.txt
+++ b/testsuite/CMakeLists.txt
@@ -19,4 +19,5 @@ Project(matmod.tests C Fortran)
 
 add_subdirectory(matmod.tests)
 add_subdirectory(Prog.tests)
+add_subdirectory(Analysis.tests)
 enable_testing()