diff --git a/Docs/source/Basic_Tutorial.rst b/Docs/source/Basic_Tutorial.rst
index 10ea65e3..6b6c058d 100644
--- a/Docs/source/Basic_Tutorial.rst
+++ b/Docs/source/Basic_Tutorial.rst
@@ -72,7 +72,8 @@ that see the Building_libamrex chapter of the amrex docs here. `here <https://am
 
 The HeatEquation examples solve a 2D or 3D (determined by how you set DIM in the GNUmakefile)
 heat equation explicitly on a domain-decomposed mesh.  This example is described in detail in
-the Basics_ chapter of the amrex Documentation
+the Basics_ chapter of the amrex Documentation.
+There is also an example code, HeadEquation_Restart, that shows how to add restart capability to a time-advancement code.
 
 .. _Basics: https://amrex-codes.github.io/amrex/docs_html/Basics.html
 
diff --git a/ExampleCodes/Basic/HeatEquation_Restart/CMakeLists.txt b/ExampleCodes/Basic/HeatEquation_Restart/CMakeLists.txt
new file mode 100644
index 00000000..6c75ade6
--- /dev/null
+++ b/ExampleCodes/Basic/HeatEquation_Restart/CMakeLists.txt
@@ -0,0 +1,15 @@
+if (AMReX_SPACEDIM EQUAL 1)
+   return()
+endif ()
+
+# List of source files
+set(_sources main.cpp checkpoint.cpp myfunc.H)
+list(TRANSFORM _sources PREPEND "Source/")
+
+# List of input files
+file( GLOB_RECURSE _input_files LIST_DIRECTORIES false Exec/input* )
+
+setup_tutorial(_sources _input_files)
+
+unset( _sources )
+unset( _input_files   )
diff --git a/ExampleCodes/Basic/HeatEquation_Restart/Exec/GNUmakefile b/ExampleCodes/Basic/HeatEquation_Restart/Exec/GNUmakefile
new file mode 100644
index 00000000..53a0d790
--- /dev/null
+++ b/ExampleCodes/Basic/HeatEquation_Restart/Exec/GNUmakefile
@@ -0,0 +1,18 @@
+# AMREX_HOME defines the directory in which we will find all the AMReX code.
+AMREX_HOME ?= ../../../../../amrex
+
+DEBUG        = FALSE
+USE_MPI      = FALSE
+USE_OMP      = FALSE
+COMP         = gnu
+DIM          = 3
+
+include $(AMREX_HOME)/Tools/GNUMake/Make.defs
+
+include ../Source/Make.package
+VPATH_LOCATIONS  += ../Source
+INCLUDE_LOCATIONS += ./Source
+
+include $(AMREX_HOME)/Src/Base/Make.package
+
+include $(AMREX_HOME)/Tools/GNUMake/Make.rules
diff --git a/ExampleCodes/Basic/HeatEquation_Restart/Exec/inputs b/ExampleCodes/Basic/HeatEquation_Restart/Exec/inputs
new file mode 100644
index 00000000..9c9e86ab
--- /dev/null
+++ b/ExampleCodes/Basic/HeatEquation_Restart/Exec/inputs
@@ -0,0 +1,11 @@
+n_cell = 32
+max_grid_size = 16
+
+nsteps = 1000
+plot_int = 100
+chk_int = 100
+
+dt = 1.e-5
+
+restart = -1
+#restart = 200
diff --git a/ExampleCodes/Basic/HeatEquation_Restart/Source/Make.package b/ExampleCodes/Basic/HeatEquation_Restart/Source/Make.package
new file mode 100644
index 00000000..9faa1837
--- /dev/null
+++ b/ExampleCodes/Basic/HeatEquation_Restart/Source/Make.package
@@ -0,0 +1,3 @@
+CEXE_sources += main.cpp
+CEXE_sources += checkpoint.cpp
+CEXE_headers += myfunc.H
diff --git a/ExampleCodes/Basic/HeatEquation_Restart/Source/checkpoint.cpp b/ExampleCodes/Basic/HeatEquation_Restart/Source/checkpoint.cpp
new file mode 100644
index 00000000..3157d488
--- /dev/null
+++ b/ExampleCodes/Basic/HeatEquation_Restart/Source/checkpoint.cpp
@@ -0,0 +1,121 @@
+#include "myfunc.H"
+
+namespace {
+    void GotoNextLine (std::istream& is)
+    {
+        constexpr std::streamsize bl_ignore_max { 100000 };
+        is.ignore(bl_ignore_max, '\n');
+    }
+}
+
+// create a checkpoint directory
+// write out time and BoxArray to a Header file
+// write out multifab data
+void WriteCheckpoint(const int& step,
+                     const Real& time,
+                     const MultiFab& phi)
+{
+
+
+    // checkpoint file name, e.g., chk00010
+    const std::string& checkpointname = Concatenate("chk",step,5);
+
+    BoxArray ba = phi.boxArray();
+
+    // single level problem
+    int nlevels = 1;
+
+    // ---- prebuild a hierarchy of directories
+    // ---- dirName is built first.  if dirName exists, it is renamed.  then build
+    // ---- dirName/subDirPrefix_0 .. dirName/subDirPrefix_nlevels-1
+    // ---- if callBarrier is true, call ParallelDescriptor::Barrier()
+    // ---- after all directories are built
+    // ---- ParallelDescriptor::IOProcessor() creates the directories
+    PreBuildDirectorHierarchy(checkpointname, "Level_", nlevels, true);
+
+    VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size);
+
+    // write Header file to store time and BoxArray
+    if (ParallelDescriptor::IOProcessor()) {
+
+        std::ofstream HeaderFile;
+        HeaderFile.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size());
+        std::string HeaderFileName(checkpointname + "/Header");
+        HeaderFile.open(HeaderFileName.c_str(), std::ofstream::out   |
+                        std::ofstream::trunc |
+                        std::ofstream::binary);
+
+        if( !HeaderFile.good()) {
+            FileOpenFailed(HeaderFileName);
+        }
+
+        HeaderFile.precision(17);
+
+        // write out title line
+        HeaderFile << "Checkpoint file for MagneX\n";
+
+        // write out time
+        HeaderFile << time << "\n";
+
+        // write the BoxArray
+        ba.writeOn(HeaderFile);
+        HeaderFile << '\n';
+    }
+
+    // write the MultiFab data to, e.g., chk00010/Level_0/
+    VisMF::Write(phi, MultiFabFileFullPrefix(0, checkpointname, "Level_", "phi"));
+}
+
+// read in the time and BoxArray, then create a DistributionMapping
+// Define phi and fill it with data from the checkpoint file
+void ReadCheckpoint(const int& restart,
+                    Real& time,
+                    MultiFab& phi,
+                    BoxArray& ba,
+                    DistributionMapping& dm)
+{
+
+    // checkpoint file name, e.g., chk00010
+    const std::string& checkpointname = Concatenate("chk",restart,5);
+
+    Print() << "Restart from checkpoint " << checkpointname << "\n";
+
+    VisMF::IO_Buffer io_buffer(VisMF::GetIOBufferSize());
+
+    std::string line, word;
+
+    // Header
+    {
+        std::string File(checkpointname + "/Header");
+        Vector<char> fileCharPtr;
+        ParallelDescriptor::ReadAndBcastFile(File, fileCharPtr);
+        std::string fileCharPtrString(fileCharPtr.dataPtr());
+        std::istringstream is(fileCharPtrString, std::istringstream::in);
+
+        // read in title line
+        std::getline(is, line);
+
+        // read in time
+        is >> time;
+        GotoNextLine(is);
+
+        // read in BoxArray from Header
+        ba.readFrom(is);
+        GotoNextLine(is);
+
+        // create a distribution mapping
+        dm.define(ba, ParallelDescriptor::NProcs());
+
+        // Nghost = number of ghost cells for each array
+        int Nghost = 1;
+
+        // Ncomp = number of components for each array
+        int Ncomp = 1;
+
+        phi.define(ba, dm, Ncomp, Nghost);
+    }
+
+    // read in the MultiFab data
+    VisMF::Read(phi, MultiFabFileFullPrefix(0, checkpointname, "Level_", "phi"));
+
+}
diff --git a/ExampleCodes/Basic/HeatEquation_Restart/Source/main.cpp b/ExampleCodes/Basic/HeatEquation_Restart/Source/main.cpp
new file mode 100644
index 00000000..bc128ba7
--- /dev/null
+++ b/ExampleCodes/Basic/HeatEquation_Restart/Source/main.cpp
@@ -0,0 +1,252 @@
+#include <AMReX_PlotFileUtil.H>
+#include <AMReX_ParmParse.H>
+
+#include "myfunc.H"
+
+using namespace amrex;
+
+int main (int argc, char* argv[])
+{
+    amrex::Initialize(argc,argv);
+
+    main_main();
+
+    amrex::Finalize();
+    return 0;
+}
+
+void main_main ()
+{
+
+    // **********************************
+    // SIMULATION PARAMETERS
+
+    // number of cells on each side of the domain
+    int n_cell;
+
+    // size of each box (or grid)
+    int max_grid_size;
+
+    // total steps in simulation
+    int nsteps;
+
+    // time step
+    Real dt;
+
+    // how often to write a plotfile
+    int plot_int;
+
+    // how often to write a checkpoint file
+    int chk_int;
+
+    // restart = -1 --> start new simulation
+    // restart >  0 --> restart from step=restart
+    int restart;
+
+    // inputs parameters
+    {
+        // ParmParse is way of reading inputs from the inputs file
+        // pp.get means we require the inputs file to have it
+        // pp.query means we optionally need the inputs file to have it - but we must supply a default here
+        ParmParse pp;
+
+        // We need to get n_cell from the inputs file - this is the number of cells on each side of
+        //   a square (or cubic) domain.
+        pp.get("n_cell",n_cell);
+
+        // The domain is broken into boxes of size max_grid_size
+        pp.get("max_grid_size",max_grid_size);
+
+        // Default nsteps to 10, allow us to set it to something else in the inputs file
+        nsteps = 10;
+        pp.query("nsteps",nsteps);
+
+        // time step
+        pp.get("dt",dt);
+
+        // Default plot_int to -1, allow us to set it to something else in the inputs file
+        //  If plot_int < 0 then no plot files will be written
+        plot_int = -1;
+        pp.query("plot_int",plot_int);
+
+        // Default chk_int to -1, allow us to set it to something else in the inputs file
+        //  If chk_int < 0 then no chk files will be written
+        chk_int = -1;
+        pp.query("chk_int",chk_int);
+
+        // restart = -1 --> start new simulation
+        // restart >  0 --> restart from step=restart
+        restart = -1;
+        pp.query("restart",restart);
+
+    }
+
+    // **********************************
+    // SIMULATION SETUP
+
+    // make Geometry object
+    // geom contains information such as the physical domain size,
+    //               number of points in the domain, and periodicity
+    Geometry geom;
+
+    // AMREX_D_DECL means "do the first X of these, where X is the dimensionality of the simulation"
+    IntVect dom_lo(AMREX_D_DECL(       0,        0,        0));
+    IntVect dom_hi(AMREX_D_DECL(n_cell-1, n_cell-1, n_cell-1));
+
+    // Make a single box that is the entire domain
+    Box domain(dom_lo, dom_hi);
+
+    // This defines the physical box, [0,1] in each direction.
+    RealBox real_box({AMREX_D_DECL( 0., 0., 0.)},
+                     {AMREX_D_DECL( 1., 1., 1.)});
+
+    // periodic in all direction
+    Array<int,AMREX_SPACEDIM> is_periodic{AMREX_D_DECL(1,1,1)};
+
+    // This defines a Geometry object
+    geom.define(domain, real_box, CoordSys::cartesian, is_periodic);
+
+    // extract dx from the geometry object
+    GpuArray<Real,AMREX_SPACEDIM> dx = geom.CellSizeArray();
+
+    // create a BoxArray and DistributionMapping
+    // ba will contain a list of boxes that cover the domain
+    // if restarting, BoxArray is read in from the checkpoint directory
+    BoxArray ba;
+    DistributionMapping dm;
+
+    // simulation time
+    Real time;
+
+    // start step - will be 1 unless restarting
+    int start_step;
+
+    MultiFab phi_old;
+    MultiFab phi_new;
+
+    // Nghost = number of ghost cells for each array
+    int Nghost = 1;
+
+    // Ncomp = number of components for each array
+    int Ncomp = 1;
+
+    // start a new simulation
+    if (restart <= 0) {
+
+        // Initialize the boxarray "ba" from the single box "domain"
+        ba.define(domain);
+
+        // Break up boxarray "ba" into chunks no larger than "max_grid_size" along a direction
+        ba.maxSize(max_grid_size);
+
+        // How Boxes are distrubuted among MPI processes
+        dm.define(ba);
+
+        // we allocate two phi multifabs; one will store the old state, the other the new.
+        phi_old.define(ba, dm, Ncomp, Nghost);
+        phi_new.define(ba, dm, Ncomp, Nghost);
+
+        // time = starting time in the simulation
+        time = 0.0;
+
+        start_step = 1;
+
+        // **********************************
+        // INITIALIZE DATA
+
+        // loop over boxes
+        for (MFIter mfi(phi_old); mfi.isValid(); ++mfi)
+        {
+            const Box& bx = mfi.validbox();
+
+            const Array4<Real>& phiOld = phi_old.array(mfi);
+
+            // set phi = 1 + e^(-(r-0.5)^2)
+            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
+            {
+                Real x = (i+0.5) * dx[0];
+                Real y = (j+0.5) * dx[1];
+#if (AMREX_SPACEDIM == 2)
+                Real rsquared = ((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5))/0.01;
+#elif (AMREX_SPACEDIM == 3)
+                Real z= (k+0.5) * dx[2];
+                Real rsquared = ((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)+(z-0.5)*(z-0.5))/0.01;
+#endif
+                phiOld(i,j,k) = 1. + std::exp(-rsquared);
+            });
+        }
+
+        // Write a plotfile of the initial data if plot_int > 0
+        if (plot_int > 0)
+        {
+            int step = 0;
+            const std::string& pltfile = amrex::Concatenate("plt",step,5);
+            WriteSingleLevelPlotfile(pltfile, phi_old, {"phi"}, geom, time, 0);
+        }
+
+    } else {
+
+        // restart a simulation
+        // read in the time and BoxArray, then create a DistributionMapping
+        // Define phi_old and fill it with data from the checkpoint file
+        ReadCheckpoint(restart, time, phi_old, ba, dm);
+
+        // set start_step to restart+1
+        start_step = restart + 1;
+
+        phi_new.define(ba, dm, Ncomp, Nghost);
+    }
+
+    for (int step = start_step; step <= nsteps; ++step)
+    {
+        // fill periodic ghost cells
+        phi_old.FillBoundary(geom.periodicity());
+
+        // new_phi = old_phi + dt * Laplacian(old_phi)
+        // loop over boxes
+        for ( MFIter mfi(phi_old); mfi.isValid(); ++mfi )
+        {
+            const Box& bx = mfi.validbox();
+
+            const Array4<Real>& phiOld = phi_old.array(mfi);
+            const Array4<Real>& phiNew = phi_new.array(mfi);
+
+            // advance the data by dt
+            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
+            {
+                phiNew(i,j,k) = phiOld(i,j,k) + dt *
+                    ( (phiOld(i+1,j,k) - 2.*phiOld(i,j,k) + phiOld(i-1,j,k)) / (dx[0]*dx[0])
+                     +(phiOld(i,j+1,k) - 2.*phiOld(i,j,k) + phiOld(i,j-1,k)) / (dx[1]*dx[1])
+#if (AMREX_SPACEDIM == 3)
+                     +(phiOld(i,j,k+1) - 2.*phiOld(i,j,k) + phiOld(i,j,k-1)) / (dx[2]*dx[2])
+#endif
+                        );
+            });
+        }
+
+        // update time
+        time = time + dt;
+
+        // copy new solution into old solution
+        MultiFab::Copy(phi_old, phi_new, 0, 0, 1, 0);
+
+        // Tell the I/O Processor to write out which step we're doing
+        amrex::Print() << "Advanced step " << step << "\n";
+
+        // Write a plotfile of the current data (plot_int was defined in the inputs file)
+        if (plot_int > 0 && step%plot_int == 0)
+        {
+            const std::string& pltfile = amrex::Concatenate("plt",step,5);
+            WriteSingleLevelPlotfile(pltfile, phi_new, {"phi"}, geom, time, step);
+        }
+
+        // Write a checkpoint file of the current data (chk_int was defined in the inputs file)
+        if (chk_int > 0 && step%chk_int == 0)
+        {
+            // create a checkpoint directory
+            // write out time and BoxArray to a Header file
+            // write out multifab data
+            WriteCheckpoint(step,time,phi_new);
+        }
+    }
+}
diff --git a/ExampleCodes/Basic/HeatEquation_Restart/Source/myfunc.H b/ExampleCodes/Basic/HeatEquation_Restart/Source/myfunc.H
new file mode 100644
index 00000000..d711199a
--- /dev/null
+++ b/ExampleCodes/Basic/HeatEquation_Restart/Source/myfunc.H
@@ -0,0 +1,20 @@
+#include "AMReX_PlotFileUtil.H"
+
+#ifndef MYFUNC_H_
+#define MYFUNC_H_
+
+using namespace amrex;
+
+void main_main ();
+
+void WriteCheckpoint(const int& step,
+                     const Real& time,
+                     const MultiFab& phi);
+
+void ReadCheckpoint(const int& restart,
+                    Real& time,
+                    MultiFab& phi,
+                    BoxArray& ba,
+                    DistributionMapping& dm);
+
+#endif