diff --git a/README.md b/README.md
index b457a6ec78943627309e0bd718a33d1496d3a4d3..c0ae292569bdc5f7f32f69d6b974fc4bb1c7ac45 100644
--- a/README.md
+++ b/README.md
@@ -1,5 +1,5 @@
 # Overview
-This project contains supplemental code for the article:
+This project contains supplemental [Julia](https://julialang.org/) code for the article:
 
 >M. Thelander Andrén, "Using Radial Basis Functions to Appriximate the LQG-Optimal Event-Based Sampling Policy",
 Submitted to European Control Conference, 2019
@@ -34,8 +34,102 @@ First make sure that you have installed Julia v0.6 and the required packages. Th
 Note that it may take quite some time to run the demo the first time, since all the required packages are then being compiled. It takes roughly 1 min to run the entire script with $`30^2=900`$ collocation points.
 
 ## Printouts
+While running the demo, the current status is printed to the Julia REPL. First, the script generates the set of collocation points:
+```Julia
+GENERATING COLLOCATION POINTS...DONE!
+```
+If specified, it will then proceed by doing a binary search for the smallest feasible value of the shape parameter $`c`$:
+```Julia
+SEARCHING FOR SMALLEST FEASIBLE C...
+cᵤ = 160.0, cₗ = 80.0
+cᵤ = 120.0, cₗ = 80.0
+cᵤ = 100.0, cₗ = 80.0
+cᵤ = 90.0, cₗ = 80.0
+cᵤ = 90.0, cₗ = 85.0
+cᵤ = 87.5, cₗ = 85.0
+cᵤ = 86.25, cₗ = 85.0
+cᵤ = 86.25, cₗ = 85.625
+cᵤ = 85.9375, cₗ = 85.625
+cᵤ = 85.9375, cₗ = 85.78125
+cᵤ = 85.859375, cₗ = 85.78125
+cᵤ = 85.8203125, cₗ = 85.78125
+cᵤ = 85.8203125, cₗ = 85.80078125
+cᵤ = 85.810546875, cₗ = 85.80078125
+FOUND C = 85.810546875
+```
+It then proceeds by generating a linear complementarity problem (LCP) to compute the RBF weights:
+```Julia
+GENERATING LCP...DONE!
+```
+With the LCP defined, it will call the QP-solver in the package `OSQP` to solve it:
+```Julia
+SOLVING LCP...
+-----------------------------------------------------------------
+           OSQP v0.4.1  -  Operator Splitting QP Solver
+              (c) Bartolomeo Stellato,  Goran Banjac
+        University of Oxford  -  Stanford University 2018
+-----------------------------------------------------------------
+problem:  variables n = 900, constraints m = 1800
+          nnz(P) + nnz(A) = 1215861
+settings: linear system solver = qdldl,
+          eps_abs = 1.0e-09, eps_rel = 1.0e-09,
+          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
+          rho = 1.00e-01 (adaptive),
+          sigma = 1.00e-06, alpha = 1.60, max_iter = 100000
+          check_termination: on (interval 25),
+          scaling: on, scaled_termination: off
+          warm start: on, polish: off
+
+iter   objective    pri res    dua res    rho        time
+   1  -2.1255e+01   1.81e+01   3.25e+00   1.00e-01   1.50e+00s
+ 200  -1.1622e+00   2.61e-01   3.72e-03   1.00e-01   2.49e+00s
+ 400  -3.0957e-01   6.63e-02   1.11e-03   1.00e-01   3.39e+00s
+ 600  -8.8080e-02   1.72e-02   2.97e-04   1.00e-01   4.34e+00s
+ 800  -2.9303e-02   4.88e-03   8.66e-05   1.00e-01   5.50e+00s
+1000  -1.0692e-02   1.55e-03   3.35e-05   1.00e-01   6.49e+00s
+1200  -4.3451e-03   5.77e-04   1.46e-05   1.00e-01   7.66e+00s
+1400  -1.7629e-03   2.28e-04   6.05e-06   1.00e-01   8.71e+00s
+1600  -7.1483e-04   9.20e-05   2.47e-06   1.00e-01   9.82e+00s
+1800  -2.8981e-04   3.72e-05   1.00e-06   1.00e-01   1.08e+01s
+2000  -1.1749e-04   1.51e-05   4.06e-07   1.00e-01   1.17e+01s
+2200  -4.7627e-05   6.12e-06   1.65e-07   1.00e-01   1.28e+01s
+2400  -1.9307e-05   2.48e-06   6.68e-08   1.00e-01   1.38e+01s
+2600  -7.8269e-06   1.00e-06   2.71e-08   1.00e-01   1.47e+01s
+2800  -3.1729e-06   4.07e-07   1.10e-08   1.00e-01   1.57e+01s
+3000  -1.2862e-06   1.65e-07   4.45e-09   1.00e-01   1.66e+01s
+3200  -5.2142e-07   6.70e-08   1.80e-09   1.00e-01   1.76e+01s
+3400  -2.1138e-07   2.71e-08   7.31e-10   1.00e-01   1.85e+01s
+3600  -8.5688e-08   1.10e-08   2.96e-10   1.00e-01   1.94e+01s
+3800  -3.4736e-08   4.46e-09   1.20e-10   1.00e-01   2.04e+01s
+3850  -2.7717e-08   3.56e-09   9.59e-11   1.00e-01   2.06e+01s
+
+status:               solved
+number of iterations: 3850
+optimal objective:    -0.0000
+run time:             2.06e+01s
+optimal rho estimate: 8.75e-02
+
+LCP SOLVED! OBJECTIVE VALUE = -2.7716763106866438e-8
+COMPUTING RBF-WEIGHTS...DONE!
+```
+Finally, it evaluates the RBF approximation on a grid for plotting. When specified
+(and $`A=0`$), it will also evaluate the true analytic solution on the same grid
+for comparison.
+```Julia
+EVALUATING RBF-APPROXIMATION ON GRID...DONE!
+EVALUATING TRUE SOLUTION ON GRID...DONE!
+MAE ON GRID = 0.0008456961574835812
+GENERATING PLOT...DONE!
+```
 
 ## Plots
+If specified, the demo script will output plots of the resulting approximation.
+Below is an example:
+![Example plot output](figures/example.png)
+
+When validation is specified, the plot will also feature a comparison to the 
+true solution:
+![Validation plot output](figures/validation.png)
 
 # Modifying the Demo Code
 We encourage you to open the file 'demo.jl' and modify it to try out different precision of the approximation and different systems. Whenever you have changed the parameters to your liking, simply run 'include('demo.jl')' in the Julia REPL again to see the new result.