diff --git a/MOT Capture Process Simulation/+Plotting/plotDynamicalQuantities.m b/MOT Capture Process Simulation/+Plotting/plotDynamicalQuantities.m
new file mode 100644
index 0000000..83f22e0
--- /dev/null
+++ b/MOT Capture Process Simulation/+Plotting/plotDynamicalQuantities.m	
@@ -0,0 +1,107 @@
+function plotDynamicalQuantities(obj, MaximumVelocity, IncidentAtomDirection, IncidentAtomPosition)
+
+    f_h = Helper.getFigureByTag('Trajectories Plot');
+    set(groot,'CurrentFigure',f_h);
+    a_h = get(f_h, 'CurrentAxes');
+    if ~isempty(get(a_h, 'Children')) 
+        clf(f_h);
+    end
+    f_h.Name  = 'Trajectories Plot';
+    f_h.Units = 'pixels';
+    
+    set(0,'units','pixels');
+    screensize = get(0,'ScreenSize');
+    f_h.Position = [[screensize(3)/50 screensize(4)/10] 1870 850];
+    
+    N     = obj.NumberOfAtoms;
+    Theta = IncidentAtomDirection;
+    z     = IncidentAtomPosition;
+    
+    obj.setInitialConditions();
+    L     = obj.OvenDistance * 2;
+    T     = obj.SimulationTime;
+    tau   = obj.TimeStep;
+    
+    Y                    = linspace(0,MaximumVelocity,N);
+    DynamicalQuantities  = zeros(length(Y),int64(T/tau),6);
+    for i=1:length(Y)
+        x =-L/2;
+        vx = Y(i)*cos(Theta);
+        vz = Y(i)*sin(Theta);
+        r  = [x,0,z];
+        v  = [vx,0,vz];
+        DynamicalQuantities(i,:,:) = obj.solver(r, v);
+    end
+
+    t = linspace(0, T, T/tau) * 1e+3;
+    for i=1:length(Y)
+        subplot(3, 3, 1)
+        hold on
+        plot(t, DynamicalQuantities(i,:,1) * 1e+3,'r','linewidth',1.3)
+        hXLabel_1 = xlabel('Time (ms)');
+        hYLabel_1 = ylabel('x (mm)');
+        hold off
+        subplot(3, 3, 2)
+        hold on
+        plot(t, DynamicalQuantities(i,:,2) * 1e+3,'g','linewidth',1.3)
+        hXLabel_2 = xlabel('Time (ms)');
+        hYLabel_2 = ylabel('y (mm)');
+        hold off
+        subplot(3, 3, 3)
+        hold on
+        plot(t, DynamicalQuantities(i,:,3) * 1e+3,'b','linewidth',1.3)
+        hXLabel_3 = xlabel('Time (ms)');
+        hYLabel_3 = ylabel('z (m/s)');
+        hold off
+        subplot(3, 3, 4)
+        hold on
+        plot(t, DynamicalQuantities(i,:,4),'r','linewidth',1.3)
+        hXLabel_4 = xlabel('Time (ms)');
+        hYLabel_4 = ylabel('v_x (m/s)');
+        hold off
+        subplot(3, 3, 5)
+        hold on
+        plot(t, DynamicalQuantities(i,:,5),'g','linewidth',1.3)
+        hXLabel_5 = xlabel('Time (ms)');
+        hYLabel_5 = ylabel('v_y (m/s)');
+        hold off
+        subplot(3, 3, 6)
+        hold on
+        plot(t, DynamicalQuantities(i,:,6),'b','linewidth',1.3)
+        hXLabel_6 = xlabel('Time (ms)');
+        hYLabel_6 = ylabel('v_z (m/s)');
+        hold off
+        subplot(3, 3, 7)
+        hold on
+        plot(DynamicalQuantities(i,:,1), DynamicalQuantities(i,:,4),'r','linewidth',1.3)
+        hXLabel_7 = xlabel('x (mm)');
+        hYLabel_7 = ylabel('v_x (m/s)');
+        xlim([-L/2 L/2])
+        hold off
+        subplot(3, 3, 8)
+        hold on
+        plot(DynamicalQuantities(i,:,2), DynamicalQuantities(i,:,5),'g','linewidth',1.3)
+        hXLabel_8 = xlabel('y (mm)');
+        hYLabel_8 = ylabel('v_y (m/s)');
+        xlim([-1 1] * 1e-04)
+        hold off
+        subplot(3, 3, 9)
+        hold on
+        plot(DynamicalQuantities(i,:,3), DynamicalQuantities(i,:,6),'b','linewidth',1.3)
+        hXLabel_9 = xlabel('z (mm)');
+        hYLabel_9 = ylabel('v_z (m/s)');
+        xlim([-1 1] * 1e-04)
+        hold off
+    end
+    
+    hold off
+    hTitle = sgtitle(sprintf("Magnetic gradient = %.2f T/m", obj.MagneticGradient));
+    set([hXLabel_1, hXLabel_2, hXLabel_3, hXLabel_4, hXLabel_5, hXLabel_6, hXLabel_7, hXLabel_8, hXLabel_9,...
+        hYLabel_1, hYLabel_2, hYLabel_3, hYLabel_4, hYLabel_5, hYLabel_6, hYLabel_7, hYLabel_8, hYLabel_9], ...
+        'FontSize'   , 14          );
+    set( hTitle                    , ...
+        'FontSize'   , 18          );
+
+    Helper.bringFiguresWithTagInForeground();
+end
+