diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 2115b297cacaee0ef7ae7010b7945e11c246a974..1807b07f45f5d02ef709f8b43d6976a1efd21689 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -365,8 +365,12 @@ int main(int argc, char *argv[]) {
     multigridStep.ignoreNodes_ = &ignoreNodes;
     // }}}
 
+    std::fstream octave_writer("data", std::fstream::out);
     timer.reset();
     auto const timesteps = parset.get<size_t>("timesteps");
+    octave_writer << "# name: A" << std::endl << "# type: matrix" << std::endl
+                  << "# rows: " << timesteps << std::endl << "# columns: 3"
+                  << std::endl;
     double const h = 1.0 / timesteps;
     for (size_t run = 1; run <= timesteps; ++run) {
       if (parset.get<bool>("printProgress")) {
@@ -468,6 +472,10 @@ int main(int argc, char *argv[]) {
                                       "%|40t|u[%03d] = %+3e");
         std::cout << boost::format(formatter) % run % (*s4_new)[i] % run % u4[i]
                   << std::endl;
+        octave_writer << (*s4_new)[i] << " " << u4[i][0] * 1e6 << " "
+                      << ((h * run <= 0.5) ? sin(h * run * 2 * M_PI) * 1
+                                           : (h * run - 0.5) * 2 * 1)
+                      << std::endl;
         break;
       }
 
@@ -552,6 +560,7 @@ int main(int argc, char *argv[]) {
       std::cout << "sup |u1 - u4| = " << diff4.infinity_norm() << ", "
                 << "|u1 - u4| = " << diff4.two_norm() << std::endl;
     }
+    octave_writer.close();
 
     if (parset.get<bool>("printFrictionalBoundary")) {
       // Print displacement on frictional boundary
diff --git a/src/plot_displacement.m b/src/plot_displacement.m
new file mode 100644
index 0000000000000000000000000000000000000000..63132b011e3614947731c047a11116e34e3d6ea2
--- /dev/null
+++ b/src/plot_displacement.m
@@ -0,0 +1,4 @@
+close all;
+
+plot(1:length(A),A(:,2));
+axis tight;
diff --git a/src/plot_neumann.m b/src/plot_neumann.m
new file mode 100644
index 0000000000000000000000000000000000000000..8fd9e217da81b4db4d93b1147fca087a36f805c9
--- /dev/null
+++ b/src/plot_neumann.m
@@ -0,0 +1,4 @@
+close all;
+
+plot(1:length(A),A(:,3));
+axis tight;
diff --git a/src/plot_state.m b/src/plot_state.m
new file mode 100644
index 0000000000000000000000000000000000000000..0172ebc79f665f369a29973b0baebc4d9acdbfc4
--- /dev/null
+++ b/src/plot_state.m
@@ -0,0 +1,4 @@
+close all;
+
+plot(1:length(A),A(:,1));
+axis tight;