diff --git a/SU2_CFD/include/solvers/CHeatSolver.hpp b/SU2_CFD/include/solvers/CHeatSolver.hpp
index a286f9be771..df2d99149b7 100644
--- a/SU2_CFD/include/solvers/CHeatSolver.hpp
+++ b/SU2_CFD/include/solvers/CHeatSolver.hpp
@@ -196,6 +196,17 @@ class CHeatSolver final : public CScalarSolver<CHeatVariable> {
                         unsigned short iMesh,
                         unsigned short iRKStep) override;
 
+ /*!
+   * \brief Source term computation.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] numerics_container - Description of the numerical method.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] iMesh - Index of the mesh in multigrid computations.
+   */
+  void Source_Residual(CGeometry *geometry, CSolver **solver_container,  CNumerics **numerics_container,
+                       CConfig *config, unsigned short iMesh) override ;
+
 
   void Set_Heatflux_Areas(CGeometry *geometry, CConfig *config) override;
 
diff --git a/SU2_CFD/src/solvers/CHeatSolver.cpp b/SU2_CFD/src/solvers/CHeatSolver.cpp
index 54769769002..479e81e42b4 100644
--- a/SU2_CFD/src/solvers/CHeatSolver.cpp
+++ b/SU2_CFD/src/solvers/CHeatSolver.cpp
@@ -311,6 +311,20 @@ void CHeatSolver::Viscous_Residual(CGeometry *geometry, CSolver **solver_contain
   }
 }
 
+void CHeatSolver::Source_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container,
+                                  CConfig *config, unsigned short iMesh) {
+
+  /*--- Regular source terms go here. ---*/
+  /*--- ... ---*/
+
+ /*--- Custom user defined source term (from the python wrapper) ---*/
+  if (config->GetPyCustomSource()) {
+    CustomSourceResidual(geometry, solver_container, numerics_container, config, iMesh);
+  }
+
+}
+
+
 void CHeatSolver::Set_Heatflux_Areas(CGeometry *geometry, CConfig *config) {
 
   BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py
index 9b8dfafe8ec..07a72635302 100755
--- a/TestCases/parallel_regression.py
+++ b/TestCases/parallel_regression.py
@@ -1534,6 +1534,17 @@ def main():
     pywrapper_zimont.command = TestCase.Command("mpirun -np 2", "python", "run.py")
     test_list.append(pywrapper_zimont)
 
+
+    # Heat solver unsteady with source 
+    pywrapper_Unst_Heat_Source = TestCase('pywrapper_Unst_Heat_Source') 
+    pywrapper_Unst_Heat_Source.cfg_dir = "py_wrapper/custom_heat_source" 
+    pywrapper_Unst_Heat_Source.cfg_file = "run.py" 
+    pywrapper_Unst_Heat_Source.test_iter = 10 
+    pywrapper_Unst_Heat_Source.unsteady = True 
+    pywrapper_Unst_Heat_Source.test_vals = [-5.235402, 300.010000, 300.000000]
+    pywrapper_Unst_Heat_Source.command = TestCase.Command("mpirun -n 2", "python", "run.py")
+    test_list.append(pywrapper_Unst_Heat_Source)
+    
     ##############################################
     ### Method of Manufactured Solutions (MMS) ###
     ##############################################
diff --git a/TestCases/py_wrapper/custom_heat_source/run.py b/TestCases/py_wrapper/custom_heat_source/run.py
new file mode 100644
index 00000000000..563278aff0c
--- /dev/null
+++ b/TestCases/py_wrapper/custom_heat_source/run.py
@@ -0,0 +1,160 @@
+#!/usr/bin/env python
+
+## \file run.py
+#  \brief Unsteady heat transfer case with custom source term.
+#  \version 8.3.0 "Harrier"
+#
+# SU2 Project Website: https://su2code.github.io
+#
+# The SU2 Project is maintained by the SU2 Foundation
+# (http://su2foundation.org)
+#
+# Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md)
+#
+# SU2 is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License, or (at your option) any later version.
+#
+# SU2 is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+
+from mpi4py import MPI
+import sys
+import pysu2			            # imports the SU2 wrapped module
+
+common_settings = """
+SOLVER= HEAT_EQUATION
+INC_NONDIM= DIMENSIONAL
+
+TIME_DOMAIN= YES
+TIME_STEP= 0.002
+TIME_MARCHING= DUAL_TIME_STEPPING-2ND_ORDER
+
+FREESTREAM_TEMPERATURE= 300
+MATERIAL_DENSITY= 1000
+SPECIFIC_HEAT_CP = 500
+THERMAL_CONDUCTIVITY_CONSTANT= 10.0
+
+MARKER_HEATFLUX= ( x_minus, 0, x_plus, 0, y_minus, 0, y_plus, 0 )
+MARKER_PYTHON_CUSTOM= ( x_minus, x_plus, y_minus, y_plus )
+MARKER_MONITORING= ( x_minus, x_plus, y_minus, y_plus )
+
+PYTHON_CUSTOM_SOURCE= YES
+
+NUM_METHOD_GRAD= GREEN_GAUSS
+CFL_NUMBER= 100
+
+LINEAR_SOLVER= CONJUGATE_GRADIENT
+DISCADJ_LIN_SOLVER= CONJUGATE_GRADIENT
+LINEAR_SOLVER_PREC= ILU
+DISCADJ_LIN_PREC= ILU
+LINEAR_SOLVER_ERROR= 1E-5
+LINEAR_SOLVER_ITER= 100
+
+MESH_FORMAT= RECTANGLE
+MESH_BOX_SIZE= ( 33, 33, 0 )
+MESH_BOX_LENGTH= ( __SIZE__, __SIZE__, 0 )
+
+OUTPUT_FILES= RESTART, PARAVIEW
+OUTPUT_WRT_FREQ= 100,100
+OBJECTIVE_FUNCTION= AVG_TEMPERATURE
+
+INNER_ITER= 20
+CONV_RESIDUAL_MINVAL= -4
+CONV_STARTITER= 2
+MAX_TIME= 10.0
+TIME_ITER= 101
+"""
+
+primal_settings = """
+MATH_PROBLEM= DIRECT
+SCREEN_OUTPUT= TIME_ITER, CUR_TIME, INNER_ITER, RMS_RES, LINSOL, AVG_TEMPERATURE, TAVG_AVG_TEMPERATURE
+HISTORY_OUTPUT= ITER, RMS_RES, HEAT, TAVG_HEAT
+"""
+
+def HeatSource(time, driver, iPoint):
+  """Applies a source in the lower left quadrant for a period of time."""
+  allCoords = driver.Coordinates()
+  coord = allCoords.Get(iPoint)
+  x = coord[0]
+  y = coord[1]
+  xp = 0.0125
+  yp = 0.0125
+  R = 0.0125
+  Source = 0.0
+  if (time > 0.0):
+      if ((x-xp)*(x-xp) + (y-yp)*(y-yp) < R*R):
+          Source = 100.0
+
+  return Source
+
+def RunPrimal(size):
+  """
+  Run the heat solver with a custom heat source.
+  Returns the final average boundary temperature.
+  """
+  comm = MPI.COMM_WORLD
+  rank = comm.Get_rank()
+
+  if rank == 0:
+    with open('config_unsteady.cfg', 'w') as f:
+      f.write(common_settings.replace('__SIZE__', str(size)) + primal_settings)
+    print("\n------------------------------ Begin Solver -----------------------------\n")
+  sys.stdout.flush()
+
+  comm.Barrier()
+
+ # Initialize the corresponding driver of SU2, this includes solver preprocessing
+  try:
+    driver = pysu2.CSinglezoneDriver("config_unsteady.cfg", 1, comm)
+  except TypeError as exception:
+    print('A TypeError occured in pysu2.CDriver : ',exception)
+    raise
+
+  # Run the time loop in python to vary the heat flux.
+  dt = driver.GetUnsteadyTimeStep()
+
+  iHEATSOLVER = driver.GetSolverIndices()['HEAT']
+  Source = driver.UserDefinedSource(iHEATSOLVER)
+  
+
+  for time_iter in range(driver.GetNumberTimeIter()):
+    # set the source term, per point
+    for i_node in range(driver.GetNumberNodes() - driver.GetNumberHaloNodes()):
+      # add source term:
+      S = HeatSource(time_iter * dt, driver, i_node)
+      Source.Set(i_node, 0, S)
+
+
+    driver.Preprocess(time_iter)
+
+    # Run one time iteration.
+    driver.Run()
+    driver.Postprocess()
+    driver.Update()
+
+    # Monitor the solver and output solution to file if required.
+    driver.Monitor(time_iter)
+    driver.Output(time_iter)
+
+  # Get the final average temperature.
+  avg_temperature = driver.GetOutputValue('AVG_TEMPERATURE')
+
+  # Finalize the solver and exit cleanly.
+  driver.Finalize()
+
+  return avg_temperature
+
+
+def main():
+
+  RunPrimal(0.1)
+
+if __name__ == '__main__':
+  main()