diff --git a/matlab/amdahl.m b/matlab/amdahl.m
new file mode 100644
index 0000000..b73fd3d
--- /dev/null
+++ b/matlab/amdahl.m
@@ -0,0 +1,115 @@
+%% Benchmark the simulation performance
+
+thread_numbers = 1:1:12;
+thread_atom_numbers = 10.^5;
+steps = 1e2*5;
+
+% Run once to force compilation.
+bench(10, 1, steps);
+
+thread_results = {};
+for atom_number=thread_atom_numbers
+    atom_results = {};
+    for thread_number=thread_numbers
+        tic
+        bench(thread_number, atom_number, steps);
+        time = toc;
+        atom_results{end+1} = struct('threads', thread_number, 'atoms', atom_number, 'time', time);
+    end
+    atom_results = cat(2,atom_results{:});
+    thread_results{end+1} = atom_results;
+end
+thread_results = cat(1,thread_results{:});
+save('amdahl.mat', 'thread_results', 'steps');
+
+f = @(A, p, x) A*(1-p+p./x);
+ft = fittype(f);
+fitResult = fit([thread_results(end,:).threads]', [thread_results(end,:).time]', ft);
+fitResult
+plot([thread_results(end,:).threads]', [thread_results(end,:).time]', 'ok'); hold on;
+xs = [thread_results(end,:).threads]';
+xs = linspace(min(xs), max(xs), 1000);
+plot(xs, fitResult(xs), 'k--'); hold off;
+
+%%
+% Plot a graph showing the results
+
+set(gcf, 'Units', 'centimeters');
+pos = get(gcf, 'Position');
+set(gcf, 'Position', [ pos(1) pos(2) 9 12 ]);
+clf;
+axes('Units', 'centimeters', 'Position', [ 1.2 4.7 7.5 7 ]);
+
+c1 = [ 0.1608 0.5804 0.6980 ];
+c0 = 0*[ 0.0118 0.0196 0.1176 ];
+get_color = @(n) interp1([0; log10(max(thread_atom_numbers))], [ c0; c1 ], log10(n));
+
+set(gcf, 'Color', 'w');
+h = [];
+for i=1:size(thread_results, 1)
+   h(i) = plot([thread_results(i,:).threads], 1e6*[thread_results(i,:).time]./([thread_results(i,:).atoms].*steps), '.-', 'Color', get_color(thread_results(i,1).atoms)); hold on;
+   %plot([thread_results(i,:).threads], [thread_results(i,:).time], '.-',
+   %'Color', get_color(thread_results(i,1).atoms)); hold on;
+   %plot([thread_results(i,:).threads], [thread_results(i,1).time]./
+   %[thread_results(i,:).threads], '--', 'Color',
+   %get_color(thread_results(i,1).atoms)); hold on;
+end
+xlabel('', 'interpreter', 'latex', 'FontSize', 11);
+% tau = total wall time per atom, per thread
+ylabel('$\tau$ ($\mu$s) ', 'Interpreter', 'latex', 'FontSize', 11);
+grid on;
+set(gca, 'GridLineStyle', ':');
+xlim([min(thread_numbers) max(thread_number)]);
+set(get(gca, 'XAxis'), 'TickLabelInterpreter', 'Latex');
+set(get(gca, 'YAxis'), 'TickLabelInterpreter', 'Latex');
+set(gca, 'YScale', 'log');
+set(gca, 'XTick', []);
+set(gca, 'YScale', 'log');
+xlim([1 6]);
+labels = arrayfun(@(x) [num2str(x) ' atoms'], thread_atom_numbers, 'UniformOutput', 0);
+selected = [ 1 5 length(labels) ];
+legend(h(selected),labels{selected}, 'Interpreter', 'Latex');
+
+% fit Amdahl's law and show on the graph
+ax2 = axes('Units', 'centimeters', 'Position', [ 1.2 1.2 7.5 3.3 ]);
+f = @(A, p, x) A*(1-p+p./x);
+ft = fittype(f);
+fitResult = fit([thread_results(end,:).threads]', [thread_results(end,:).time]', ft);
+plot([thread_results(end,:).threads]', [thread_results(end,:).time]', 'o', 'Color', c1); hold on;
+xs = [thread_results(end,:).threads]';
+xs = linspace(min(xs), max(xs), 1000);
+plot(xs, fitResult(xs), 'k--'); hold off;
+grid on;
+set(ax2, 'GridLineStyle', ':');
+set(gca, 'XTick', 1:6);
+set(get(ax2, 'XAxis'), 'TickLabelInterpreter', 'Latex');
+set(get(ax2, 'YAxis'), 'TickLabelInterpreter', 'Latex');
+xlabel('number of threads', 'interpreter', 'latex', 'FontSize', 11);
+% tau = total wall time per atom, per thread
+ylabel('wall time (s) ', 'Interpreter', 'latex', 'FontSize', 11);
+
+% Render to file
+%set(gcf, 'Units', 'centimeters');
+%pos = get(gcf, 'Position');
+%set(gcf, 'Position', [ pos(1) pos(2) 9 7.5 ]);
+pos = get(gcf, 'Position');
+w = pos(3); 
+h = pos(4);
+p = 0.01;
+set(gcf,...
+  'PaperUnits','centimeters',...
+  'PaperPosition',[p*w p*h w h],...
+  'PaperSize',[w*(1+2*p) h*(1+2*p)]);
+set(gcf, 'Renderer', 'painters')
+saveas(gcf, 'bench.pdf')
+
+function bench(thread_number, atom_numbers, steps)
+   
+config = struct('n_threads', int32(thread_number), 'n_steps',  int32(steps), 'n_atoms',  int32(atom_numbers));
+fH = fopen('benchmark.json', 'w');
+fprintf(fH, '%s', jsonencode(config));
+fclose(fH);
+
+system('cargo run --example benchmark --release');
+
+end
\ No newline at end of file
diff --git a/src/ecs.rs b/src/ecs.rs
index 8658a22..8e7a872 100644
--- a/src/ecs.rs
+++ b/src/ecs.rs
@@ -15,7 +15,6 @@ use crate::integrator::{EulerIntegrationSystem, Step};
 use crate::laser;
 use crate::laser::repump::Dark;
 use crate::magnetic;
-use crate::optimization::LargerEarlyTimestepOptimizationSystem;
 use crate::output::console_output::ConsoleOutputSystem;
 use crate::sim_region;
 use specs::{Dispatcher, DispatcherBuilder, World};
@@ -26,7 +25,6 @@ pub fn register_components(world: &mut World) {
 	magnetic::register_components(world);
 	laser::register_components(world);
 	atom_sources::register_components(world);
-	//detector::register_components(world);
 	sim_region::register_components(world);
 	world.register::<Dark>();
 }
@@ -54,7 +52,6 @@ pub fn create_simulation_dispatcher_builder() -> DispatcherBuilder<'static, 'sta
 			"add_gravity",
 		],
 	);
-	//builder = detector::add_systems_to_dispatch(builder, &[]);
 	builder = builder.with(ConsoleOutputSystem, "", &["euler_integrator"]);
 	builder = builder.with(DeleteToBeDestroyedEntitiesSystem, "", &["euler_integrator"]);
 	builder = sim_region::add_systems_to_dispatch(builder, &[]);
diff --git a/src/laser/photons_scattered.rs b/src/laser/photons_scattered.rs
index 00e7f20..259b68d 100644
--- a/src/laser/photons_scattered.rs
+++ b/src/laser/photons_scattered.rs
@@ -259,13 +259,13 @@ impl<'a> System<'a> for CalculateActualPhotonsScatteredSystem {
 
         match fluctuations_option {
             None => {
-                for (expected, actual) in
-                    (&expected_photons_vector, &mut actual_photons_vector).join()
-                {
-                    for index in 0..expected.contents.len() {
-                        actual.contents[index].scattered = expected.contents[index].scattered;
-                    }
-                }
+                (&expected_photons_vector, &mut actual_photons_vector)
+                    .par_join()
+                    .for_each(|(expected, actual)| {
+                        for index in 0..expected.contents.len() {
+                            actual.contents[index].scattered = expected.contents[index].scattered;
+                        }
+                    });
             }
             Some(_rand) => {
                 (&expected_photons_vector, &mut actual_photons_vector)
diff --git a/src/laser/repump.rs b/src/laser/repump.rs
index 4f3dd18..21cdb48 100644
--- a/src/laser/repump.rs
+++ b/src/laser/repump.rs
@@ -39,14 +39,17 @@ impl<'a> System<'a> for RepumpSystem {
         Entities<'a>,
     );
     fn run(&mut self, (repump_opt, lazy, num, ent): Self::SystemData) {
+        use rayon::prelude::*;
+        use specs::ParJoin;
+
         match repump_opt {
             None => (),
             Some(repump) => {
-                for (ent, num) in (&ent, &num).join() {
+                (&ent, &num).par_join().for_each(|(ent, num)| {
                     if repump.if_loss(num.total) {
                         lazy.insert(ent, Dark {})
                     }
-                }
+                });
             }
         }
     }