diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m
index e06c41f..f09a03e 100644
--- a/Classes/00_signals/Signal.m
+++ b/Classes/00_signals/Signal.m
@@ -172,11 +172,9 @@ classdef Signal
 
             hold on;
             if isempty(options.color)
-                % plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1);
-                plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1);
+                plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1);

             else
-                % plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1,'Color',options.color);
-                plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Color',options.color);
+                plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1,'Color',options.color);

             end
             % 2 c)
             % - xlabel if not already here: time in readable format (1 ms and not 1e-3 s)
diff --git a/Classes/02_optical/DP_Fiber.m b/Classes/02_optical/DP_Fiber.m
index a5f3dce..c1f08ea 100644
--- a/Classes/02_optical/DP_Fiber.m
+++ b/Classes/02_optical/DP_Fiber.m
@@ -24,6 +24,8 @@ classdef DP_Fiber
         SS_dzmax        % [m]   max dz (adaptive SSFM)
         SS_dzmin        % [m]   min dz (adaptive SSFM)
         n_waveplates    % number of PMD waveplates
+        useGPU          % GPU acceleration: true, false, or 'auto' (default)
+        useSingle       % Use single precision on GPU (default: false)
 
         % ---- Internal state (persistent between calls) ----
         state           % struct mirroring legacy 'state'
@@ -56,6 +58,8 @@ classdef DP_Fiber
                 options.SS_dzmax     = 2e4                     % m
                 options.SS_dzmin     = 100                     % m
                 options.n_waveplates = 100
+                options.useGPU       = 'auto'          % 'auto', true, or false
+                options.useSingle    = false           % single precision GPU
             end
 
             % Copy provided options into properties
@@ -208,7 +212,7 @@ classdef DP_Fiber
                 % Frequency-dependent PMD phase term (legacy form)
                 st.brf.db0     = (R.rand(st.wave_plates,1)*2*pi - pi) * brf_multiplier;
                 st.brf.db1     = sqrt(3*pi/8)*(st.dgd/obj.fa)/st.wave_plates .* st.omega;
-                st.brf.simdgd  = 0; 
+                st.brf.simdgd  = 0;
                 % cumsum used in legacy only for debug; keep compatibility variable:
                 ~cumsum(st.brf.db0); % no-op to mirror legacy path
 
@@ -228,7 +232,18 @@ classdef DP_Fiber
             x_in = signal_in(:,1).';
             y_in = signal_in(:,2).';
 
-            [x_out, y_out, obj.state] = CNLSE_plain(x_in, y_in, obj.state);
+            % Determine GPU usage
+            if ischar(obj.useGPU) || isstring(obj.useGPU)
+                if strcmpi(obj.useGPU, 'auto')
+                    gpuFlag = [];  % Let CNLSE_plain auto-detect
+                else
+                    error('DP_Fiber:InvalidGPU', 'useGPU must be true, false, or ''auto''');
+                end
+            else
+                gpuFlag = logical(obj.useGPU);
+            end
+
+            [x_out, y_out, obj.state] = CNLSE_plain(x_in, y_in, obj.state, gpuFlag, obj.useSingle);
 
             obj.state.propagated_length = obj.state.propagated_length + obj.state.L;
 
diff --git a/Classes/02_optical/Optical_Demultiplex.m b/Classes/02_optical/Optical_Demultiplex.m
index 7d7c39e..8e272d2 100644
--- a/Classes/02_optical/Optical_Demultiplex.m
+++ b/Classes/02_optical/Optical_Demultiplex.m
@@ -42,7 +42,7 @@ classdef Optical_Demultiplex < handle
 
         function signalclasses_out = process(obj, signalclass_in)
 
-            % ---- Infer wavelength: either given or from input total signal 
+            % ---- Infer wavelength: either given or from input total signal
             if isempty(obj.wavelengthplan)
                 obj.wavelengthplan = signalclass_in.lambda; %meter
             else
@@ -81,7 +81,7 @@ classdef Optical_Demultiplex < handle
                 obj
                 signal_in
             end
-            
+
             w = obj.fs_out ./ obj.fs_in ;
             blocklen_in = length(signal_in);
             blocklen_out = w*blocklen_in;
@@ -119,30 +119,31 @@ classdef Optical_Demultiplex < handle
             N = size(lo,1);
             C = size(lo,2);
 
-            x_envelopes = zeros(N, C, 'like', signal_in);
-            y_envelopes = zeros(N, C, 'like', signal_in);
+            % ---- VECTORIZED: Process all channels in parallel ----
+            % Batched FFT operates on each column simultaneously on GPU
 
-            s1 = signal_in(:,1);
-            s2 = signal_in(:,2);
+            % Extract polarization signals
+            s1 = signal_in(:,1);  % X polarization [N×1]
+            s2 = signal_in(:,2);  % Y polarization [N×1]
 
-            % Reusable work buffers (avoid reallocations)
-            wrk_time = zeros(N,1, 'like', signal_in);
-            wrk_freq = zeros(N,1, 'like', signal_in);
+            % Broadcast signal to all channels and multiply with LO
+            % s1, s2 are [N×1], lo is [N×C] → result is [N×C]
+            x_mixed = att .* s1 .* lo;  % [N×C]
+            y_mixed = att .* s2 .* lo;  % [N×C]
+
+            % Batched FFT: each column computed in parallel
+            x_freq = fft(x_mixed);  % [N×C]
+            y_freq = fft(y_mixed);  % [N×C]
+
+            % Apply filter (H is [N×1], broadcasts across columns)
+            x_filtered = x_freq .* H;  % [N×C]
+            y_filtered = y_freq .* H;  % [N×C]
+
+            % Batched IFFT
+            x_envelopes = ifft(x_filtered);  % [N×C]
+            y_envelopes = ifft(y_filtered);  % [N×C]
 
-            for c = 1:C
-                % ---- X branch ----
-                wrk_time(:) = att .* s1 .* lo(:,c);   % N×1
-                wrk_freq(:) = fft(wrk_time);          % N×1
-                wrk_freq(:) = wrk_freq .* H;          % N×1
-                x_envelopes(:,c) = ifft(wrk_freq);     % N×1
 
-                % ---- Y branch ----
-                wrk_time(:) = att .* s2 .* lo(:,c);
-                wrk_freq(:) = fft(wrk_time);
-                wrk_freq(:) = wrk_freq .* H;
-                y_envelopes(:,c) = ifft(wrk_freq);
-            end
- 
         end
     end
 end
diff --git a/Classes/02_optical/Optical_Multiplex.m b/Classes/02_optical/Optical_Multiplex.m
index 1d722b0..7fc8a33 100644
--- a/Classes/02_optical/Optical_Multiplex.m
+++ b/Classes/02_optical/Optical_Multiplex.m
@@ -1,10 +1,10 @@
 classdef Optical_Multiplex < handle
     % Ta