Plotting the time-depending function in the parabolic case: 
********************************************************************

Open a file by  

  string filename="datafield.dat";
  ofstream file;
  file.open(filename.c_str());  

Then, you can use  

  void printTimeStep(double ti, const Vector& xnodes, const Vector& u_ti, ostream& file){
    for(int j=0;j<xnodes.size();j++)
      file<<ti<<" "<<xnodes[j]<<" "<<u_ti[j]<<endl;
  }

(where ti is the current time, xnodes contains the nodes of the space-mesh, and u_ti is the 
solution vector at the current time step) to write each data point at time ti to the file 
whenever you computed the solution vector at a time step.

At the end, don't forget to close the file again by 
  
  file.close();

Then, the file should have the correct format to serve as input for Mathematica:
Create a Mathematica notebook in the same folder as your data file, and do something like 
  
  SetDirectory[NotebookDirectory[]];
  data = Import["datafield.dat"];
  plot = ListPlot3D[data (*add some options, as eg. AxesLabel, PlotRange,..*)];
  Export["datafield.jpg", plot];
  Print[plot];


Alternatively, you can go for the corresponding function in matlab.