diff --git a/examples/2d_hertz_contact.jl b/examples/2d_hertz_contact.jl
index 4dcbe7f..8ad1ad8 100644
--- a/examples/2d_hertz_contact.jl
+++ b/examples/2d_hertz_contact.jl
@@ -24,11 +24,6 @@
 # ``a = 6.21 \;\mathrm{mm}``.
 
 using JuliaFEM
-using JuliaFEM.Preprocess
-using JuliaFEM.Postprocess
-using Logging
-Logging.configure(level=INFO)
-add_elements! = JuliaFEM.add_elements!
 
 # Simulation starts by reading the mesh. Model is constructed and meshed using
 # SALOME, thus mesh format is .med. Mesh type is quite simple structure,
@@ -135,18 +130,15 @@ add_slave_elements!(contact, contact_slave_elements)
 # results using ParaView, thus we write our results to Xdmf format, which uses
 # well defined standards XML and HDF to store model data.
 
-step = Analysis(Nonlinear)
-add_problems!(step, [upper, lower, bc_fixed, bc_sym_23, load, contact])
+analysis = Analysis(Nonlinear)
+add_problems!(analysis, upper, lower, bc_fixed, bc_sym_23, load, contact)
 xdmf = Xdmf("2d_hertz_results"; overwrite=true)
-## todo for 2d
-## for body in (upper, lower)
-##     push!(body.postprocess_fields, "stress")
-## end
-add_results_writer!(step, xdmf)
+add_results_writer!(analysis, xdmf)
 
 # In last part, we run the analysis.
 
-step()
+run!(analysis)
+close(xdmf)
 
 # # Results
 
@@ -172,7 +164,7 @@ end
 
 println("2d hertz contact resultant forces: Rn = $Rn, Rt = $Rt")
 
-using Base.Test
+using Test
 @test isapprox(Rn, 35.0e3)
 @test isapprox(Rt, 0.0)
 
@@ -217,5 +209,3 @@ end
 println("Contact radius: $a_rad")
 
 # This example briefly described some of the core features of JuliaFEM.
-
-close(xdmf.hdf) # src
diff --git a/examples/3d_frame.jl b/examples/3d_frame.jl
index 58f6abb..8bcfb38 100644
--- a/examples/3d_frame.jl
+++ b/examples/3d_frame.jl
@@ -11,12 +11,7 @@
 
 # ![](3d_frame/model.png)
 
-using JuliaFEM
-using JuliaFEM.Preprocess
-using FEMBase.Test
-using Logging
-Logging.configure(level=INFO)
-add_elements! = JuliaFEM.add_elements!
+using JuliaFEM, LinearAlgebra
 
 # Reading mesh
 
@@ -32,7 +27,7 @@ println("Number of nodes in a model: ", length(mesh.nodes))
 # system and polar moment of inertia.
 
 beam_elements = create_elements(mesh, "FRAME")
-info("Number of elements: ", length(beam_elements))
+@info("Number of elements: ", length(beam_elements))
 update!(beam_elements, "youngs modulus", 210.0e6)
 update!(beam_elements, "shear modulus", 84.0e6)
 update!(beam_elements, "density", 7850.0e-3)
@@ -50,7 +45,7 @@ update!(beam_elements, "polar moment of inertia", 30.0e-5)
 for element in beam_elements
     X1, X2 = element("geometry", 0.0)
     t = (X2-X1)/norm(X2-X1)
-    I = eye(3)
+    I = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
     k = indmax([norm(cross(t, I[:,k])) for k in 1:3])
     n = cross(t, I[:,k])/norm(cross(t, I[:,k]))
     update!(element, "normal", n)
@@ -77,12 +72,12 @@ add_elements!(frame, bc_elements)
 
 # Perform modal analysis
 
-step = Analysis(Modal)
+analysis = Analysis(Modal)
 xdmf = Xdmf(joinpath(datadir, "3d_frame_results"); overwrite=true)
-add_results_writer!(step, xdmf)
-add_problems!(step, [frame])
-run!(step)
-close(xdmf.hdf)
+add_results_writer!(analysis, xdmf)
+add_problems!(analysis, frame)
+run!(analysis)
+close(xdmf)
 
 # Each `Analysis` can have properties, e.g. time, maximum number of iterations,
 # convergence tolerance and so on. Eigenvalues of calculation are stored as a
diff --git a/examples/generate_stiffness_matrices.jl b/examples/generate_stiffness_matrices.jl
index 65497e9..dc9e9c4 100644
--- a/examples/generate_stiffness_matrices.jl
+++ b/examples/generate_stiffness_matrices.jl
@@ -16,4 +16,4 @@ problem = Problem(Elasticity, "test problem", 2)
 problem.properties.formulation = :plane_stress
 add_elements!(problem, [element])
 assemble!(problem, 0.0)
-K = round.(full(problem.assembly.K), 5)
+K = round.(Matrix(problem.assembly.K), 5)
diff --git a/examples/linear_static.jl b/examples/linear_static.jl
index e9bbc4a..639e372 100755
--- a/examples/linear_static.jl
+++ b/examples/linear_static.jl
@@ -6,10 +6,8 @@
 # ![](linear_static/freecad.png)
 
 # ## Preprocessing
+
 using JuliaFEM
-using JuliaFEM.Preprocess
-using JuliaFEM.Postprocess
-add_elements! = JuliaFEM.add_elements!
 
 # First we will read in the mesh. Geometry and mesh are greated with FreeCAD,
 # where med format is selected for exporting. Mesh file consist also edge and
@@ -46,7 +44,8 @@ filter!(is_not_Seg3_or_Tri6, model.elements)
 # define a function, which finds nodes on the given plane yz, xz or xy from the
 # given height.
 
-function add_nodes_at_certain_plane_to_node_set!(mesh, name, vector_id, distance, radius=6.0)
+function add_nodes_at_certain_plane_to_node_set!(mesh, name, vector_id, distance,
+                                                 radius=6.0)
     for (node, coords) in mesh.nodes
         if isapprox(coords[vector_id], distance, atol=radius)
             add_node_to_node_set!(mesh, name, node)
@@ -58,8 +57,10 @@ end
 # We will find nodes from the xz-plane going through point (0,50,0) or actually
 # we previously defined the radius to be 6.0, which means (0,[44,56],0). In other
 # words we will select each node, which second coordinate value is between 44
-# and 56. This function will edit mesh object and add node set called `:mid_fixed` to it. 
-add_nodes_at_certain_plane_to_node_set!(mesh,:mid_fixed,2,50.0)
+# and 56. This function will edit mesh object and add node set called `:mid_fixed`
+# to it. 
+
+add_nodes_at_certain_plane_to_node_set!(mesh, :mid_fixed, 2, 50.0)
 
 # We need to somehow handle the i's dot. I looked the rough coordinates of the
 # dot in FreeCAD and now we can search three closest nodes to these coordinates.
@@ -91,26 +92,19 @@ analysis = Analysis(Linear, model, fixed)
 xdmf = Xdmf("model_results"; overwrite=true)
 add_results_writer!(analysis, xdmf)
 
+# This is how the stresses are requested
+push!(model.postprocess_fields, "stress")
+
 # Now we have all we need to run the analysis.
 
 run!(analysis)
 
 # ## Postprocessing
 
-# This is how the stresses are requested
-push!(model.postprocess_fields, "stress")
-
-# Finally let's write the results to the xdmf (this is bug, the solver should
-# do this automatically, see issue #203)
-
-time = 0.0
-JuliaFEM.postprocess!(analysis, time)
-JuliaFEM.write_results!(analysis, time)
-
 # In order to look the results, we will need to close the xdmf that it is actually
 # written to the file from buffer.
 
-close(xdmf.hdf)
+close(xdmf)
 
 # Finally when we open the model in ParaView and set some settings we have this
 # end result.