diff --git a/test/simulation/sim_er.jl b/test/simulation/sim_er.jl
index 09f7698309578450ccdbd4fd9b5565a2a364c0e0..002b32a06123c12731a2a2501351631fc68e7fa9 100644
--- a/test/simulation/sim_er.jl
+++ b/test/simulation/sim_er.jl
@@ -1,14 +1,12 @@
 
 using MarkovProcesses 
-using PyPlot
+using Plots
 
 load_model("ER")
 
 σ = simulate(ER)
-plt.figure()
-plt.step(times(σ), σ[:P], "ro--", marker="x", where="post", linewidth=1.0)
-plt.savefig(get_module_path() * "/test/simulation/res_pics/sim_er.png", dpi=480)
-plt.close()
+plot(times(σ), σ[:P], linetype=:steppost, linewidth=1.0)
+savefig(get_module_path() * "/test/simulation/res_pics/sim_er.png", dpi=480)
 
 return true
 
diff --git a/test/simulation/sim_er_row_buffer_bounded.jl b/test/simulation/sim_er_row_buffer_bounded.jl
index b4c3e022a7e60ba7ba3070127343c790fb52d65a..77c490cf6650c286ebcd8be15928ba392cde1617 100644
--- a/test/simulation/sim_er_row_buffer_bounded.jl
+++ b/test/simulation/sim_er_row_buffer_bounded.jl
@@ -1,16 +1,14 @@
 
 using MarkovProcesses 
 include(get_module_path() * "/src/_tests_simulate.jl")
-using PyPlot
+using Plots
 
 load_model("_bench_perf_test/ER_row_buffer")
 ER_row_buffer.time_bound = 10.0
 
 σ = _simulate_row_buffer(ER_row_buffer)
-plt.figure()
-plt.step(times(σ), _get_values_row(σ,:P), "ro--", marker="x", where="post", linewidth=1.0)
-plt.savefig(get_module_path() * "/test/simulation/res_pics/sim_er_row_buffer_bounded.svg")
-plt.close()
+plot(times(σ), _get_values_row(σ,:P), linetype=:steppost, linewidth=1.0)
+savefig(get_module_path() * "/test/simulation/res_pics/sim_er_row_buffer_bounded.svg")
 
 return true
 
diff --git a/test/simulation/sim_sir.jl b/test/simulation/sim_sir.jl
index c0e8da5ff98a7a76a4b6f61fc144cec758f1a9da..86187dc3291738c45a55ceddbe416d07d352318e 100644
--- a/test/simulation/sim_sir.jl
+++ b/test/simulation/sim_sir.jl
@@ -1,14 +1,12 @@
 
 using MarkovProcesses 
-using PyPlot
+using Plots
 
 load_model("SIR")
 
 σ = simulate(SIR)
-plt.figure()
-plt.step(times(σ), σ[:I], "ro--", marker="x", where="post", linewidth=1.0)
-plt.savefig(get_module_path() * "/test/simulation/res_pics/sim_sir.svg")
-plt.close()
+plot(times(σ), σ[:I], linetype=:steppost, linewidth=1.0)
+savefig(get_module_path() * "/test/simulation/res_pics/sim_sir.svg")
 
 return true
 
diff --git a/test/simulation/sim_sir_bounded.jl b/test/simulation/sim_sir_bounded.jl
index 1e9ab1ff5d569734bdcdb4a21c35e84d070542a3..861177c39d7bd28f7ea52b65752afa12b44a3610 100644
--- a/test/simulation/sim_sir_bounded.jl
+++ b/test/simulation/sim_sir_bounded.jl
@@ -1,15 +1,13 @@
 
 using MarkovProcesses 
-using PyPlot
+using Plots
 
 load_model("SIR")
 SIR.time_bound = 100.0
 
 σ = simulate(SIR)
-plt.figure()
-plt.step(times(σ), σ[:I], "ro--", marker="x", where="post", linewidth=1.0)
-plt.savefig(get_module_path() * "/test/simulation/res_pics/sim_sir_bounded.svg")
-plt.close()
+plot(times(σ), σ[:I], linetype=:steppost, linewidth=1.0)
+savefig(get_module_path() * "/test/simulation/res_pics/sim_sir_bounded.svg")
 
 return true
 
diff --git a/test/simulation/sim_sir_col_buffer_bounded.jl b/test/simulation/sim_sir_col_buffer_bounded.jl
index d3f6887a52618b8e07bece852aa2d8d3e3cc8a01..5b2ccc8fc90226ed6f260392763d1bc255860f9f 100644
--- a/test/simulation/sim_sir_col_buffer_bounded.jl
+++ b/test/simulation/sim_sir_col_buffer_bounded.jl
@@ -1,16 +1,14 @@
 
 using MarkovProcesses 
 include(get_module_path() * "/src/_tests_simulate.jl")
-using PyPlot
+using Plots
 
 load_model("_bench_perf_test/SIR_col_buffer")
 SIR_col_buffer.time_bound = 100.0
 
 σ = _simulate_col_buffer(SIR_col_buffer)
-plt.figure()
-plt.step(times(σ), _get_values_col(σ,:I), "ro--", marker="x", where="post", linewidth=1.0)
-plt.savefig(get_module_path() * "/test/simulation/res_pics/sim_sir_col_buffer_bounded.svg")
-plt.close()
+plot(times(σ), _get_values_row(σ,:I), linetype=:steppost, linewidth=1.0)
+savefig(get_module_path() * "/test/simulation/res_pics/sim_sir_col_buffer_bounded.svg")
 
 return true
 
diff --git a/test/simulation/sim_sir_row_buffer_bounded.jl b/test/simulation/sim_sir_row_buffer_bounded.jl
index 329aea537f5557d955a73f9ca202b76db6c1b132..201d4f62548c8afc7f45f2d779dac756119b868f 100644
--- a/test/simulation/sim_sir_row_buffer_bounded.jl
+++ b/test/simulation/sim_sir_row_buffer_bounded.jl
@@ -1,16 +1,14 @@
 
 using MarkovProcesses 
 include(get_module_path() * "/src/_tests_simulate.jl")
-using PyPlot
+using Plots
 
 load_model("_bench_perf_test/SIR_row_buffer")
 SIR_row_buffer.time_bound = 100.0
 
 σ = _simulate_row_buffer(SIR_row_buffer)
-plt.figure()
-plt.step(times(σ), _get_values_row(σ,:I), "ro--", marker="x", where="post", linewidth=1.0)
-plt.savefig(get_module_path() * "/test/simulation/res_pics/sim_sir_row_buffer_bounded.svg")
-plt.close()
+plot(times(σ), _get_values_row(σ,:I), linetype=:steppost, linewidth=1.0)
+savefig(get_module_path() * "/test/simulation/res_pics/sim_sir_row_buffer_bounded.svg")
 
 return true