diff --git a/src/TestProblem/Hydro/Bondi/ExtAcc_Bondi.cpp b/src/TestProblem/Hydro/Bondi/ExtAcc_Bondi.cpp
index 44af678bbb..8f8ab7b011 100644
--- a/src/TestProblem/Hydro/Bondi/ExtAcc_Bondi.cpp
+++ b/src/TestProblem/Hydro/Bondi/ExtAcc_Bondi.cpp
@@ -21,6 +21,11 @@
 
 extern double Bondi_MassBH;
 extern double Bondi_Soften_R;
+extern bool   Bondi_Soliton;
+extern double Bondi_Soliton_m22;
+extern double Bondi_Soliton_rc;
+extern int    Bondi_Soliton_type;
+extern double Bondi_Soliton_t;
 
 //-------------------------------------------------------------------------------------------------------
 // Function    :  SetExtAccAuxArray_Bondi
@@ -43,6 +48,38 @@ void SetExtAccAuxArray_Bondi( double AuxArray[], const double Time )
    AuxArray[2] = amr->BoxCenter[2];
    AuxArray[3] = NEWTON_G*Bondi_MassBH;   // gravitational_constant*point_source_mass
    AuxArray[4] = Bondi_Soften_R;          // soften_length (<=0.0 --> disable)
+   AuxArray[5] = Bondi_Soliton;
+   if( Bondi_Soliton )
+   {
+      AuxArray[6] = Bondi_Soliton_m22;
+      AuxArray[7] = Bondi_Soliton_rc;
+      AuxArray[8] = NEWTON_G*Const_Msun/UNIT_M;
+      AuxArray[9] = UNIT_L/Const_kpc;
+      if( Bondi_Soliton_type > 0 )
+      {
+          switch(Bondi_Soliton_type)
+          {
+              case 1: // arctan function
+                  AuxArray[8] *= 2/(real)3.14159265*ATAN(Time/Bondi_Soliton_t);
+                  break;
+              case 2: // linear function
+                  if( Time < Bondi_Soliton_t )
+                      AuxArray[8] *= Time/Bondi_Soliton_t;
+                  break;
+              case 3: // smooth step function
+                  if( Time < Bondi_Soliton_t )
+                      AuxArray[8] *= 3*SQR(Time/Bondi_Soliton_t)-2*CUBE(Time/Bondi_Soliton_t);
+                  break;
+              case 4: // sigmoid
+                  AuxArray[8] *= 2/(1+exp(-Time*log(3)/Bondi_Soliton_t))-1;
+                  break;
+              case 5: // tanh
+                  AuxArray[8] *= tanh(Time/Bondi_Soliton_t);
+                  break;
+                
+          }
+      } 
+   }
 
 } // FUNCTION : SetExtAccAuxArray_Bondi
 #endif // #ifndef __CUDACC__
@@ -73,13 +110,31 @@ static void ExtAcc_Bondi( real Acc[], const double x, const double y, const doub
 {
 
    const double Cen[3] = { UserArray[0], UserArray[1], UserArray[2] };
-   const real GM       = (real)UserArray[3];
+         real GM       = (real)UserArray[3];
    const real eps      = (real)UserArray[4];
+   const bool SOL      = (real)UserArray[5];
    const real dx       = (real)(x - Cen[0]);
    const real dy       = (real)(y - Cen[1]);
    const real dz       = (real)(z - Cen[2]);
    const real r        = SQRT( dx*dx + dy*dy + dz*dz );
 
+   if( SOL )
+   {
+      const real m22      = (real)UserArray[6];
+      const real rc       = (real)UserArray[7];  // In code unit or kpc
+      const real Coeff    = (real)UserArray[8];
+      const real UNIT_L   = (real)UserArray[9];
+
+#ifdef Plummer
+      double M = GM*CUBE(r)/pow(SQR(r)+SQR(rc),1.5);
+#else
+      real a = SQRT(POW(2.0,1.0/8.0)-1)*(r/rc);
+      real M = (real)4.17e9/(SQR(m22/1e-1)*(rc*UNIT_L*1e3)*POW(SQR(a)+1, 7.0))*((real)3465*POW(a,13.0)+(real)23100*POW(a,11.0)+(real)65373*POW(a,9.0)+(real)101376*POW(a,7.0)+(real)92323*POW(a,5.0)+(real)48580*POW(a,3.0)-(real)3465*a+(real)3465*POW(SQR(a)+1, 7.0)*ATAN(a));
+      M *= Coeff;
+#endif
+      GM += M;
+   }
+
 // Plummer
 #  if   ( defined SOFTEN_PLUMMER )
    const real _r3 = ( eps <= (real)0.0 ) ? (real)1.0/CUBE(r) : POW( SQR(r)+SQR(eps), (real)-1.5 );
diff --git a/src/TestProblem/Hydro/Bondi/Flu_ResetByUser_Bondi.cpp b/src/TestProblem/Hydro/Bondi/Flu_ResetByUser_Bondi.cpp
index 9f8fffb30f..13a7c42ead 100644
--- a/src/TestProblem/Hydro/Bondi/Flu_ResetByUser_Bondi.cpp
+++ b/src/TestProblem/Hydro/Bondi/Flu_ResetByUser_Bondi.cpp
@@ -2,8 +2,7 @@
 
 #if ( MODEL == HYDRO  &&  defined GRAVITY )
 
-
-
+extern double Bondi_MassBH;
 extern double Bondi_InBC_Rho;
 extern double Bondi_InBC_R;
 extern double Bondi_InBC_E;
@@ -19,8 +18,8 @@ extern double Bondi_SinkEk;
 extern double Bondi_SinkEt;
 extern int    Bondi_SinkNCell;
 
-
-
+extern bool   Bondi_void;
+extern bool   Bondi_dynBH;
 
 //-------------------------------------------------------------------------------------------------------
 // Function    :  Flu_ResetByUser_Func_Bondi
@@ -52,6 +51,9 @@ int Flu_ResetByUser_Func_Bondi( real fluid[], const double Emag, const double x,
                                 const double dt, const int lv, double AuxArray[] )
 {
 
+   if ( !Bondi_void ) 
+      return false;
+
    const double Pos[3]  = { x, y, z };
    const double InBC_R2 = SQR( Bondi_InBC_R );
    double dr2[3], r2;
@@ -113,6 +115,9 @@ void Flu_ResetByUser_API_Bondi( const int lv, const int FluSg, const int MagSg,
    int    Reset;
    real   fluid[NCOMP_TOTAL], fluid_bk[NCOMP_TOTAL];
    double x, y, z, x0, y0, z0;
+// Define varible to record sink mass at every time step
+   double SinkMass_OneSubStep_ThisRank = 0;
+   double SinkMass_OneSubStep_AllRank;
 
 // reset to 0 since we only want to record the number of void cells **for one sub-step**
    Bondi_SinkNCell = 0;
@@ -120,7 +125,7 @@ void Flu_ResetByUser_API_Bondi( const int lv, const int FluSg, const int MagSg,
 
 #  pragma omp parallel for private( Reset, fluid, fluid_bk, x, y, z, x0, y0, z0 ) schedule( runtime ) \
    reduction(+:Bondi_SinkMass, Bondi_SinkMomX, Bondi_SinkMomY, Bondi_SinkMomZ, Bondi_SinkMomXAbs, Bondi_SinkMomYAbs, Bondi_SinkMomZAbs, \
-               Bondi_SinkEk, Bondi_SinkEt, Bondi_SinkNCell)
+               Bondi_SinkEk, Bondi_SinkEt, Bondi_SinkNCell, SinkMass_OneSubStep_ThisRank )
    for (int PID=0; PID<amr->NPatchComma[lv][1]; PID++)
    {
       x0 = amr->patch[0][lv][PID]->EdgeL[0] + 0.5*dh;
@@ -193,12 +198,23 @@ void Flu_ResetByUser_API_Bondi( const int lv, const int FluSg, const int MagSg,
                Bondi_SinkEk      += dv*Ek;
                Bondi_SinkEt      += dv*Et;
                Bondi_SinkNCell   ++;
+
+               SinkMass_OneSubStep_ThisRank += dv*fluid_bk[DENS];
             }
-         } // if ( Reset )
 
+            else if ( amr->patch[0][lv][PID]->son == -1 )
+            {
+//             void region must be completely refined to the max level
+               Aux_Error( ERROR_INFO, "void region lies outside the max-level region !!\n" );
+	    }
+         } // if ( Reset )
       }}} // i,j,k
    } // for (int PID=0; PID<amr->NPatchComma[lv][1]; PID++)
 
+   MPI_Allreduce( &SinkMass_OneSubStep_ThisRank, &SinkMass_OneSubStep_AllRank, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
+
+   if( Bondi_dynBH ){ Bondi_MassBH += SinkMass_OneSubStep_AllRank; }
+
 } // FUNCTION : Flu_ResetByUser_API_Bondi
 
 
diff --git a/src/TestProblem/Hydro/Bondi/Init_TestProb_Hydro_Bondi.cpp b/src/TestProblem/Hydro/Bondi/Init_TestProb_Hydro_Bondi.cpp
index f57f4af155..5175391fee 100644
--- a/src/TestProblem/Hydro/Bondi/Init_TestProb_Hydro_Bondi.cpp
+++ b/src/TestProblem/Hydro/Bondi/Init_TestProb_Hydro_Bondi.cpp
@@ -1,5 +1,8 @@
 #include "GAMER.h"
 
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_odeiv2.h>
 
 
 // problem-specific global variables
@@ -74,11 +77,32 @@ static double Bondi_HSE_Beta_Rcore;       // core radius (input parameter)
 static double Bondi_HSE_Beta_Rho0;        // peak density (set by Bondi_HSE_Dens_NormR/D)
 static double Bondi_HSE_Beta_P1;          // P(r) = P1*( 1/x + atan(x) ) + P2 assuming beta=2/3, where x=r/Rcore,
 static double Bondi_HSE_Beta_P2;          // P1=G*MassBH*Rho0/Rcore, and P2 currently fixed to -0.5*pi*P1 so that P(inf)=0
+
+// parameters for soliton
+       bool   Bondi_void;
+       bool   Bondi_dynBH;
+       bool   Bondi_Soliton;
+       double Bondi_Soliton_m22;
+       int    Bondi_Soliton_type;
+       double Bondi_Soliton_t;
+       double Bondi_Soliton_rc;
+       double Bondi_Soliton_MassHalo;
+       double Bondi_Soliton_Redshift;
+static double *Bondi_Soliton_PresProf[2] = { NULL, NULL };   // pressure profile tabke: [0/1] = [radius/density]
+
+// parameters for bondi initial condition
+       bool   Bondi_Init;
+static char   Bondi_Init_Filename[1000];
+static double *Bondi_Init_Data=NULL;
+static int    Bondi_Init_bin;
+       double *Bondi_Init_Prof_r, *Bondi_Init_Prof_d, *Bondi_Init_Prof_v, *Bondi_Init_Prof_p;
+
 // =======================================================================================
 
 
 // problem-specific function prototypes
 void Init_ExtAcc_Bondi();
+void Init_ExtPot_Bondi();
 void Record_Bondi();
 bool Flag_Bondi( const int i, const int j, const int k, const int lv, const int PID, const double *Threshold );
 int Flu_ResetByUser_Func_Bondi( real fluid[], const double Emag, const double x, const double y, const double z, const double Time,
@@ -89,6 +113,34 @@ static void BondiBC( real Array[], const int ArraySize[], real fluid[], const in
                      const int GhostSize, const int idx[], const double pos[], const double Time,
                      const int lv, const int TFluVarIdxList[], double AuxArray[] );
 
+void SetExtAccAuxArray_Bondi( double [], const double );
+
+
+
+
+#ifdef GRAVITY
+//-------------------------------------------------------------------------------------------------------
+// Function    :  Poi_UserWorkBeforePoisson_Bondi
+// Description :  Call SetExtAccAuxArray_Bondi() to reset Bondi_MassBH before invoking the Poisson solver
+//
+// Note        :  1. Invoked by Gra_AdvanceDt() using the function pointer "Poi_UserWorkBeforePoisson_Ptr"
+//
+// Parameter   :  Time : Target physical time
+//                lv   : Target refinement level
+//
+// Return      :  None
+//////-------------------------------------------------------------------------------------------------------
+void Poi_UserWorkBeforePoisson_Bondi( const double Time, const int lv )
+{
+
+   SetExtAccAuxArray_Bondi( ExtAcc_AuxArray, Time );
+
+#  ifdef GPU
+   CUAPI_SetConstMemory_ExtAccPot();
+#  endif
+
+} // FUNCTION : Poi_UserWorkBeforePoisson_Bondi
+#endif // #ifdef GRAVITY
 
 
 
@@ -212,6 +264,20 @@ void SetParameter()
    ReadPara->Add( "Bondi_HSE_Pres_NormT", &Bondi_HSE_Pres_NormT,      false,        Useless_bool,     Useless_bool      );
    ReadPara->Add( "Bondi_HSE_Beta_Rcore", &Bondi_HSE_Beta_Rcore,     -1.0,          Eps_double,       NoMax_double      );
 
+   ReadPara->Add( "Bondi_void",               &Bondi_void,                    true,         Useless_bool,     Useless_bool      );
+   ReadPara->Add( "Bondi_dynBH",              &Bondi_dynBH,                   true,         Useless_bool,     Useless_bool      );
+   ReadPara->Add( "Bondi_Soliton",            &Bondi_Soliton,                 false,        Useless_bool,     Useless_bool      );
+   ReadPara->Add( "Bondi_Soliton_m22",        &Bondi_Soliton_m22,            -1.0,          NoMin_double,     NoMax_double      );
+   ReadPara->Add( "Bondi_Soliton_type",       &Bondi_Soliton_type,            0,          NoMin_int,        NoMax_int         );
+   ReadPara->Add( "Bondi_Soliton_t",          &Bondi_Soliton_t,              -1.0,          NoMin_double,     NoMax_double      );
+   ReadPara->Add( "Bondi_Soliton_rc",         &Bondi_Soliton_rc,             -1.0,          NoMin_double,     NoMax_double      );
+   ReadPara->Add( "Bondi_Soliton_MassHalo",   &Bondi_Soliton_MassHalo,       -1.0,          NoMin_double,     NoMax_double      );
+   ReadPara->Add( "Bondi_Soliton_Redshift",   &Bondi_Soliton_Redshift,       -1.0,          NoMin_double,     NoMax_double      );
+
+
+   ReadPara->Add( "Bondi_Init",               &Bondi_Init,                    false,        Useless_bool,     Useless_bool      );
+   ReadPara->Add( "Bondi_Init_Filename",      Bondi_Init_Filename,           Useless_str,  Useless_str,      Useless_str       );
+
    ReadPara->Read( FileName );
 
    delete ReadPara;
@@ -328,6 +394,50 @@ void SetParameter()
       }
    } // if ( Bondi_HSE )
 
+   if ( Bondi_Soliton )
+   {
+      if( Bondi_Soliton_rc < 0.0 ){
+         double z = Bondi_Soliton_Redshift;
+         double Mh = Bondi_Soliton_MassHalo;
+         double H0 = 67.66;
+         double Om0 = 0.3111;
+         H0 = H0*1e5/(Const_kpc*1e3);
+         double a0 = Om0*SQR(H0)/(2.47e-5*SQR(1e7/(Const_kpc*1e3)));
+         double H_H0_z = 1/(1-Om0*(1-pow(1+z,3.0)*(1+z+a0)/a0));
+         double Om_z   = Om0*pow(1+z,3.0)*H_H0_z;
+         double Om_0   = Om0/(1-Om0*(1-(1+a0)/a0));
+         double kiz_z  = (18*SQR(3.14159265)+82*(Om_z-1)-39*SQR(Om_z-1))/Om_z;
+         double kiz_0  = (18*SQR(3.14159265)+82*(Om_0-1)-39*SQR(Om_0-1))/Om_0;
+         
+         Bondi_Soliton_rc = 1.6/Bondi_Soliton_m22*pow(1+z,-1.0/2.0)*pow(kiz_z/kiz_0,-1.0/6.0)*pow(Mh/1e9,-1.0/3.0);
+      }
+      Bondi_Soliton_rc *= UnitExt_L/UNIT_L;
+      Bondi_Soliton_t  *= Bondi_TimeB;
+      Soliton_SetPresProfileTable();
+   }
+
+   if ( Bondi_Init )
+   {
+            bool RowMajor_No  = false;
+      const bool AllocMem_Yes = true;
+      const int  NCol         = 4;// Radius, Density, Velocity, MachNumber
+      const int  Col[NCol]    = {1,3,5,6};
+      Bondi_Init_bin = Aux_LoadTable( Bondi_Init_Data, Bondi_Init_Filename, NCol, Col, RowMajor_No, AllocMem_Yes );
+
+      // Convert to code unit and mach number->pressure
+      Bondi_Init_Prof_r = Bondi_Init_Data + 0*Bondi_Init_bin;
+      Bondi_Init_Prof_d = Bondi_Init_Data + 1*Bondi_Init_bin;
+      Bondi_Init_Prof_v = Bondi_Init_Data + 2*Bondi_Init_bin;
+      Bondi_Init_Prof_p = Bondi_Init_Data + 3*Bondi_Init_bin;
+      
+      for( int b=0; b<Bondi_Init_bin; b++ )
+      {
+         Bondi_Init_Prof_r[b] *= Const_kpc/UNIT_L;
+         Bondi_Init_Prof_d[b] *= 5.0e-25; // in cgs
+         Bondi_Init_Prof_v[b] *= 1000*1e5;  // in cm/s
+         Bondi_Init_Prof_p[b] = SQR(Bondi_Init_Prof_v[b]/Bondi_Init_Prof_p[b])*Bondi_Init_Prof_d[b]/GAMMA;
+      }
+    }
 
 // (4) reset other general-purpose parameters
 //     --> a helper macro PRINT_RESET_PARA is defined in Macro.h
@@ -383,11 +493,28 @@ void SetParameter()
       Aux_Message( stdout, "  Bondi_HSE_Beta        = %13.7e\n",                 Bondi_HSE_Beta                                                );
       Aux_Message( stdout, "  Bondi_HSE_Beta_Rho0   = %13.7e (%13.7e g/cm^3)\n", Bondi_HSE_Beta_Rho0, Bondi_HSE_Beta_Rho0*UNIT_D               );
       Aux_Message( stdout, "  Bondi_HSE_Beta_Rcore  = %13.7e (%13.7e kpc)\n",    Bondi_HSE_Beta_Rcore, Bondi_HSE_Beta_Rcore*UNIT_L/Const_kpc   ); }
+
+      Aux_Message( stdout, "  Bondi_void                = %s\n",                     (Bondi_void)?"YES":"NO"                                        );
+      Aux_Message( stdout, "  Bondi_dynBH               = %s\n",                     (Bondi_dynBH)?"YES":"NO"                                        );
+      Aux_Message( stdout, "  Bondi_Soliton             = %s\n",                     (Bondi_Soliton)?"YES":"NO"                                        );
+      if( Bondi_Soliton ) {
+      Aux_Message( stdout, "  Bondi_Soliton_m22         = %13.7e\n",                 Bondi_Soliton_m22                                                 );
+      Aux_Message( stdout, "  Bondi_Soliton_type        = %d\n",                     Bondi_Soliton_type                                                );
+      Aux_Message( stdout, "  Bondi_Soliton_t           = %13.7e (%13.7e Myr)\n",    Bondi_Soliton_t, Bondi_Soliton_t*UNIT_T/Const_Myr                 );
+      Aux_Message( stdout, "  Bondi_Soliton_rc          = %13.7e (%13.7e kpc)\n",    Bondi_Soliton_rc, Bondi_Soliton_rc*UNIT_L/Const_kpc               );
+      Aux_Message( stdout, "  Bondi_Soliton_MassHalo    = %13.7e Msun\n",            Bondi_Soliton_MassHalo                                            );
+      Aux_Message( stdout, "  Bondi_Soliton_Redshift    = %13.7e\n",                 Bondi_Soliton_Redshift                                            );}
+
+
+      Aux_Message( stdout, "  Bondi_Init            = %s\n",                     (Bondi_Init)?"YES":"NO"                                       );
+      if( Bondi_Init ){
+      Aux_Message( stdout, "  Bondi_Init_Filename   = %s\n",                     Bondi_Init_Filename                                           );}
+
       Aux_Message( stdout, "=============================================================================\n" );
    } // if ( MPI_Rank == 0 )
 
 
-   if ( MPI_Rank == 0 )    Aux_Message( stdout, "   Setting runtime parameters ... done\n" );
+  if ( MPI_Rank == 0 )    Aux_Message( stdout, "   Setting runtime parameters ... done\n" );
 
 } // FUNCTION : SetParameter
 
@@ -462,6 +589,41 @@ void SetGridIC( real fluid[], const double x, const double y, const double z, co
       else
          Aux_Error( ERROR_INFO, "unsupported Bondi_HSE_Mode (%d) !!\n", Bondi_HSE_Mode );
    } // if ( Bondi_HSE )
+   else if ( Bondi_Init )
+   {
+      // Set Initial Condition
+      const double r = sqrt( SQR(x-amr->BoxCenter[0]) + SQR(y-amr->BoxCenter[1]) + SQR(z-amr->BoxCenter[2]) );
+      Dens = Mis_InterpolateFromTable( Bondi_Init_bin, Bondi_Init_Prof_r, Bondi_Init_Prof_d, r )/UNIT_D;
+      Pres = Mis_InterpolateFromTable( Bondi_Init_bin, Bondi_Init_Prof_r, Bondi_Init_Prof_p, r )/UNIT_P;
+      
+      const double vmag = Mis_InterpolateFromTable( Bondi_Init_bin, Bondi_Init_Prof_r, Bondi_Init_Prof_v, r )/UNIT_V;
+      //const double dh   = amr->dh[lv];
+      //const real   dv   = CUBE(dh);
+      const double dx = x-amr->BoxCenter[0];
+      const double dy = y-amr->BoxCenter[1];
+      const double dz = z-amr->BoxCenter[2];
+      MomX = -Dens*vmag*dx/r;
+      MomY = -Dens*vmag*dy/r;
+      MomZ = -Dens*vmag*dz/r;
+      //Aux_Message( stdout, "%13.7e %13.7e %13.7e\n", MomX, MomY, MomZ );
+   }
+
+   else if ( Bondi_Init )
+   {
+      const double r = sqrt( SQR(x-amr->BoxCenter[0]) + SQR(y-amr->BoxCenter[1]) + SQR(z-amr->BoxCenter[2]) );
+      Dens = 1.945*pow( Bondi_Soliton_m22*1e1, -2.0 )*pow( Bondi_Soliton_rc*UNIT_L/Const_pc, -4.0 )*1e12/pow( 1+(9.1e-2)*SQR(r/Bondi_Soliton_rc), 8.0 );
+      Dens *= Const_Msun/CUBE(Const_pc);
+      Dens *= 1/UNIT_D;
+
+      const double *Table_R = Bondi_Soliton_PresProf[0];
+      const double *Table_P = Bondi_Soliton_PresProf[1];
+      Pres = Mis_InterpolateFromTable( 100000, Table_R, Table_P, r );
+      if( Pres==NULL_REAL ){
+         Aux_Error( ERROR_INFO, "%.3e, %.3e, %.3e\n",Table_R[0],Table_R[1000-1],r);
+         Aux_Error( ERROR_INFO, "Wrong Table\n");
+      }
+      Pres *= 1/(UNIT_P);
+   } 
 
 
 // uniform background
@@ -585,6 +747,61 @@ void HSE_SetDensProfileTable()
 } // FUNCTION : HSE_SetDensProfileTable
 
 
+//-------------------------------------------------------------------------------------------------------
+// Function    :  Soliton_SetPresProfileTable
+// Description :  Set up the pressure profile table for Soliton
+//
+// Note        :  1. Assume P(r->inf)=0
+//
+// Parameter   :  None
+//-------------------------------------------------------------------------------------------------------
+int odefunc ( double x, const double y[], double f[], void *params)
+{
+   x /= Const_kpc;
+   double rc = Bondi_Soliton_rc*UNIT_L/Const_kpc;
+   double m22 = Bondi_Soliton_m22;
+   
+   double rho = 1.945*pow(m22/1e-1, -2.0)*pow(rc*1e3, -4.0)*1e12/pow(1+(9.1e-2)*SQR(x/rc), 8.0)*Const_Msun/pow(Const_pc, 3.0);
+   double a = sqrt(pow(2.0,1.0/8.0)-1)*(x/rc);
+   double M = 4.17e9/(SQR(m22/1e-1)*(rc*1e3)*pow(SQR(a)+1, 7.0))*(3465*pow(a,13.0)+23100*pow(a,11.0)+65373*pow(a,9.0)+101376*pow(a,7.0)+92323*pow(a,5.0)+48580*pow(a,3.0)-3465*a+3465*pow(SQR(a)+1, 7.0)*atan(a))*Const_Msun;
+   f[0] = -Const_NewtonG*M*rho/SQR(x*Const_kpc);
+
+   return GSL_SUCCESS;
+}
+int * jac;
+void Soliton_SetPresProfileTable()
+{
+
+// allocate table --> deallocated by End_Bondi()
+   const int    NBin = 100000;
+   const double r_min = 0.1*amr->dh[MAX_LEVEL]*UNIT_L;
+   const double r_max = (0.5*sqrt(3.0)*amr->BoxSize[0])*UNIT_L;
+   for (int v=0; v<2; v++)    Bondi_Soliton_PresProf[v] = new double [NBin];
+   
+   int dim = 1;
+   gsl_odeiv2_system sys = {odefunc, NULL, dim, NULL};
+
+   gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk4, 1e-6, 1e-6, 0.0);
+
+   double x0 = -r_max, xf = -r_min;
+   double x = x0;
+   double y[1] = { 0.0 };
+   for( int b=1; b<=NBin; b++)
+   {
+      double xi = x0 + b*(xf-x0)/NBin;
+      int status = gsl_odeiv2_driver_apply (d, &x, xi, y);
+      if (status != GSL_SUCCESS)
+      {
+            Aux_Error( ERROR_INFO, "Error in Soliton_SetPresProfileTable, return value=%d\n", status );
+            break;
+      }
+      Bondi_Soliton_PresProf[0][NBin-b] = -x/UNIT_L;
+      Bondi_Soliton_PresProf[1][NBin-b] = y[0];
+   }
+   gsl_odeiv2_driver_free (d);
+} // void Soliton_SetDensProfileTable()
+  
+
 //-------------------------------------------------------------------------------------------------------
 // Function    :  End_Bondi
 // Description :  Free memory before terminating the program
@@ -600,6 +817,8 @@ void End_Bondi()
    {
       delete [] Bondi_HSE_DensProf[v];
       Bondi_HSE_DensProf[v] = NULL;
+      delete [] Bondi_Soliton_PresProf[v];
+      Bondi_Soliton_PresProf[v] = NULL;
    }
 
 } // FUNCTION : End_Bondi
@@ -673,7 +892,8 @@ void Init_TestProb_Hydro_Bondi()
    Flu_ResetByUser_API_Ptr  = Flu_ResetByUser_API_Bondi;
    End_User_Ptr             = End_Bondi;
 #  ifdef GRAVITY
-   Init_ExtAcc_Ptr         = Init_ExtAcc_Bondi;
+   Init_ExtAcc_Ptr               = Init_ExtAcc_Bondi;
+   Poi_UserWorkBeforePoisson_Ptr = Poi_UserWorkBeforePoisson_Bondi;
 #  endif
 #  endif // #if ( MODEL == HYDRO  &&  defined GRAVITY )
 
diff --git a/src/TestProblem/Hydro/Bondi/Record_Bondi.cpp b/src/TestProblem/Hydro/Bondi/Record_Bondi.cpp
index 8dd04bc738..5795548805 100644
--- a/src/TestProblem/Hydro/Bondi/Record_Bondi.cpp
+++ b/src/TestProblem/Hydro/Bondi/Record_Bondi.cpp
@@ -6,6 +6,7 @@
 
 // specific global variables declared in Init_TestProb_Hydro_Bondi.cpp
 // =======================================================================================
+extern double Bondi_MassBH;
 extern double Bondi_SinkMass;
 extern double Bondi_SinkMomX;
 extern double Bondi_SinkMomY;
@@ -49,10 +50,10 @@ void Record_Bondi()
          if ( Aux_CheckFileExist(FileName) )    Aux_Message( stderr, "WARNING : file \"%s\" already exists !!\n", FileName );
 
          FILE *File_User = fopen( FileName, "a" );
-         fprintf( File_User, "#%9s%16s%20s%20s%20s%20s%20s%20s%20s%20s%20s%20s%20s%20s\n",
+         fprintf( File_User, "#%9s%16s%20s%20s%20s%20s%20s%20s%20s%20s%20s%20s%20s%20s%20s\n",
                   "Step", "Time [yr]", "NVoidCell", "Mass [Msun]", "Time [yr]", "dM/dt [Msun/yr]",
                   "MomX [g*cm/s]", "MomY [g*cm/s]", "MomZ [g*cm/s]", "MomXAbs [g*cm/s]", "MomYAbs [g*cm/s]", "MomZAbs [g*cm/s]",
-                  "Ek [erg]", "Et [erg]" );
+                  "Ek [erg]", "Et [erg]", "BHMass [Msun]" );
          fclose( File_User );
       }
 
@@ -69,7 +70,7 @@ void Record_Bondi()
 
 
 //    get the total amount of sunk variables
-      double Mass_Sum, MomX_Sum, MomY_Sum, MomZ_Sum, MomXAbs_Sum, MomYAbs_Sum, MomZAbs_Sum, Ek_Sum, Et_Sum;
+      double Mass_Sum, MomX_Sum, MomY_Sum, MomZ_Sum, MomXAbs_Sum, MomYAbs_Sum, MomZAbs_Sum, Ek_Sum, Et_Sum, Mass_BH;
 
       MPI_Reduce( &Bondi_SinkMass,    &Mass_Sum,    1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
       MPI_Reduce( &Bondi_SinkMomX,    &MomX_Sum,    1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
@@ -90,6 +91,7 @@ void Record_Bondi()
       MomZAbs_Sum *= UNIT_M*UNIT_L/UNIT_T;
       Ek_Sum      *= UNIT_E;
       Et_Sum      *= UNIT_E;
+      Mass_BH      = Bondi_MassBH*UNIT_M/Const_Msun;
 
       if ( MPI_Rank == 0 )
       {
@@ -98,9 +100,9 @@ void Record_Bondi()
          dTime *= UNIT_T/Const_yr;
 
          FILE *File_User = fopen( FileName, "a" );
-         fprintf( File_User, "%10ld%16.7e%20d%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e\n",
+         fprintf( File_User, "%10ld%16.7e%20d%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e\n",
                   Step, Time[0]*UNIT_T/Const_yr, SinkNCell_Sum, Mass_Sum, dTime, Mass_Sum/dTime,
-                  MomX_Sum, MomY_Sum, MomZ_Sum, MomXAbs_Sum, MomYAbs_Sum, MomZAbs_Sum, Ek_Sum, Et_Sum );
+                  MomX_Sum, MomY_Sum, MomZ_Sum, MomXAbs_Sum, MomYAbs_Sum, MomZAbs_Sum, Ek_Sum, Et_Sum, Mass_BH );
          fclose( File_User );
       }
    } // if ( FirstTime ) ... else ...