diff --git a/classes/DataWrapper.m b/classes/DataWrapper.m
new file mode 100644
index 0000000000000000000000000000000000000000..04c827361cc842c7630607379fb4069804c775f4
--- /dev/null
+++ b/classes/DataWrapper.m
@@ -0,0 +1,94 @@
+classdef DataWrapper < handle
+    % OutputWrapper: A fast way to exchange data by pointer
+    %
+    %   Constructor:
+    %       - OutputWrapper: Creates an instance of the Wrapper class.
+    %
+    %   WARNING: Intended for data exchange where output arguments are inaccessible.
+
+    properties(SetAccess = private)
+        data        struct
+    end
+
+    properties(Access = private)
+        cache
+        counter     int32 = 0
+        dummyStrc   struct
+        chunkSize   int32                                % Size of preallocated memory
+        dataLen     int32                                % Current size of cache
+    end
+
+    methods
+        function obj = DataWrapper(dataDummy, chunkSize)
+            arguments
+                dataDummy struct
+                chunkSize {mustBeInteger, mustBePositive} = 5
+            end
+            obj.chunkSize = chunkSize;
+
+            fields = fieldnames(dataDummy)'; fields{2, 1} = [];
+            obj.dummyStrc = struct(fields{:});
+
+            obj.cache = obj.dummyStrc;
+            obj.data = obj.dummyStrc;                    % Set struct dataType
+            obj.data(obj.chunkSize) = obj.dummyStrc;     % Preallocate struct
+
+            obj.dataLen = obj.chunkSize;
+        end
+
+        function setCache(obj, data)
+            obj.cache = data;
+        end
+
+        function setData(obj)
+            obj.counter = obj.counter + 1;
+            if obj.counter > obj.dataLen
+                len = obj.dataLen + obj.chunkSize;
+                obj.dataLen = len;
+                obj.data(len) = obj.dummyStrc;
+            end
+            obj.data(obj.counter) = obj.cache;
+        end
+
+        function out = popData(obj)
+            out = obj.packageData(obj.data(1:obj.counter));
+            obj.reset();
+        end
+
+        function out = getData(obj)
+            out = obj.data(1:obj.counter);
+        end
+
+        function obj = reset(obj)
+            obj.cache = obj.dummyStrc;
+            obj.data = obj.dummyStrc;                    % Set struct dataType
+            obj.data(obj.chunkSize) = obj.dummyStrc;     % Preallocate struct
+
+            obj.dataLen = obj.chunkSize;
+            obj.counter = 0;
+        end
+    end
+
+    methods(Static)
+        function out = packageData(data)
+            % Converts struct array to struct with array fields
+            %
+            % Main assumption: fields to concatenate are colums vectors or
+            % scalars
+
+            fields = fieldnames(data);
+
+            for i = 1:length(fields)
+                currentField = fields{i};
+                currentValue = [data.(currentField)];
+
+                if isstruct(currentValue)
+                    out.(currentField) = DataWrapper.packageData([data.(currentField)]);
+                else
+                    out.(currentField) = currentValue;
+                end
+            end
+        end
+    end
+end
+
diff --git a/functions/eventFunctions/eventApogee.m b/functions/eventFunctions/eventApogee.m
index a8a8116b25bc75cbcc690e561a0af9857f954b09..e586a95b9e2de4c19e9ed4e0fcfbfca1f556eb94 100644
--- a/functions/eventFunctions/eventApogee.m
+++ b/functions/eventFunctions/eventApogee.m
@@ -1,4 +1,4 @@
-function [value, isterminal, direction] = eventApogee(t, Y, rocket, ~, ~, ~)
+function [value, isterminal, direction] = eventApogee(t, Y, rocket, varargin)
 %{
 eventApogee - Event function to stop simulation at apogee checking when a 
               value tends to zero; the value taken is to account is the 
diff --git a/functions/miscellaneous/recallOdeFcn.m b/functions/miscellaneous/recallOdeFcn.m
index 1a69a51ccf819ed220d4ef76e55f853e461e22a8..6af4a4e27a2bef4082dbb0174f8d4c23f0b2db77 100644
--- a/functions/miscellaneous/recallOdeFcn.m
+++ b/functions/miscellaneous/recallOdeFcn.m
@@ -1,100 +1,39 @@
-function [allSteps] = recallOdeFcn(fun, T, Y, varargin)
-%{
-recallOdeFcn - This function allows to compute some parameters used
-               inside the ODE integrations.
-
-INPUTS:
-        - fun, function, ode function used;
-        - T, double [n° variation, 1], integration time vector;
-        - Y, double [n° variation, 16], state matrix,
-                            [ x y z | u v w | p q r | q0 q1 q2 q3 | m | Ixx Iyy Izz ]:
-                            * (x y z), NED{north, east, down} horizontal frame;
-                            * (u v w), body frame velocities;
-                            * (p q r), body frame angular rates;
-                            * m , total mass;
-                            * (Ixx Iyy Izz), Inertias;
-                            * (q0 q1 q2 q3), attitude unit quaternion.
-
-OUTPUTS:
-        - allSteps, struct, which contains all the parameters needed.
-
-CALLED FUNCTIONS: -
-
-REVISIONS:
-- #0    16/11/2018, Release, Adriano Filippo Inno
-Copyright © 2021, Skyward Experimental Rocketry, AFD department
-All rights reserved
-
-SPDX-License-Identifier: GPL-3.0-or-later
-
-%}
-[~,firstStep] = fun(T(1),Y(:,1),varargin{:});
-
-namesFields = fieldnames(firstStep);
-NT = length(T);
-
-steps = cell(NT,1);
-for k = 1:NT
-    [~,steps{k}] = fun(T(k),Y(:,k),varargin{:});
-end
-
-for i = 1:numel(namesFields)
-    currentFieldName = namesFields{i};
-    currentField = firstStep.(currentFieldName);
-    
-    if isstruct(currentField)
-        namesArrays = fieldnames(currentField);
-        
-        for j = 1:numel(namesArrays)
-            currentArrayName = namesArrays{j};
-            currentArray = currentField.(currentArrayName);
-            sizeArray = size(currentArray);
-            sizeArrayCut = sizeArray;
-            
-            if any(sizeArray == 1)
-                sizeArrayCut = max(sizeArray);
-            end
-            
-            allSteps.(currentFieldName).(currentArrayName) = zeros([sizeArrayCut,NT]);
-            
-            for k = 1:NT
-                currentStep = steps{k};
-                
-                if all(sizeArray ~= 1)
-                    allSteps.(currentFieldName).(currentArrayName)(:,:,k) =...
-                        currentStep.(currentFieldName).(currentArrayName);
-                else
-                    allSteps.(currentFieldName).(currentArrayName)(:,k) =...
-                        currentStep.(currentFieldName).(currentArrayName);
-                end
-            end
-            
-        end
-        
-    else
-        sizeArray = size(currentField);
-        
-        sizeArrayCut = sizeArray;
-        if any(sizeArray == 1)
-            sizeArrayCut = max(sizeArray);
-        end
-        
-        allSteps.(currentFieldName) = zeros([sizeArrayCut,NT]);
-        
-        for k = 1:NT
-            currentStep = steps{k};
-            
-            if all(sizeArray ~= 1)
-                allSteps.(currentFieldName)(:,:,k) =...
-                    currentStep.(currentFieldName);
-            else
-                allSteps.(currentFieldName)(:,k) =...
-                    currentStep.(currentFieldName);
-            end
-        end
-        
-    end
-    
-end
-end
-
+function out = recallOdeFcn(fun, t, Y, varargin)
+% recallOdeFcn - This function allows to compute some parameters used
+%                inside the ODE integrations.
+% 
+% INPUTS:
+%         - fun, function, ode function used;
+%         - T, double [n° variation, 1], integration time vector;
+%         - Y, double [n° variation, 16], state matrix,
+%                             [ x y z | u v w | p q r | q0 q1 q2 q3 | m | Ixx Iyy Izz ]:
+%                             * (x y z), NED{north, east, down} horizontal frame;
+%                             * (u v w), body frame velocities;
+%                             * (p q r), body frame angular rates;
+%                             * m , total mass;
+%                             * (Ixx Iyy Izz), Inertias;
+%                             * (q0 q1 q2 q3), attitude unit quaternion.
+% 
+% OUTPUTS:
+%         - allSteps, struct, which contains all the parameters needed.
+% 
+% CALLED FUNCTIONS: -
+% 
+% REVISIONS:
+% - #0    16/11/2018, Release, Adriano Filippo Inno
+% Copyright © 2021, Skyward Experimental Rocketry, AFD department
+% All rights reserved
+% 
+% SPDX-License-Identifier: GPL-3.0-or-later
+
+wrapper = varargin{end};
+dataRaw = wrapper.getData();
+
+%% Adding t = 0, Fixing t(end)
+
+[~, data0] = fun(t(1), Y(:,1), varargin{1:end-1});
+data = [data0, dataRaw];
+
+[~, data(end)] = fun(t(end), Y(:,end), varargin{1:end-1});
+out = wrapper.packageData(data);
+wrapper.reset();
\ No newline at end of file
diff --git a/functions/odeFunctions/ballistic.m b/functions/odeFunctions/ballistic.m
index c8e56f14618858419eb488d6f9e861ad9a5ce565..c4ddb8d5cb4372596c41038e0bffa0e12b819eb7 100644
--- a/functions/odeFunctions/ballistic.m
+++ b/functions/odeFunctions/ballistic.m
@@ -1,4 +1,4 @@
-function [dY, parout] = ballistic(t, Y, rocket, environment, wind, angleRef, timeChangeRef, t0)
+function [dY, parout] = ballistic(t, Y, rocket, environment, wind, angleRef, timeChangeRef, t0, wrapper)
 arguments
     t
     Y
@@ -8,6 +8,7 @@ arguments
     angleRef        double = []
     timeChangeRef   double = []
     t0              double = 0
+    wrapper                = [] %DataWrapper = DataWrapper.empty
 end
 % ballistic - ode function of the 6DOF Rigid Rocket Model
 % 
@@ -319,31 +320,29 @@ dY = dY';
 
 %% SAVING THE QUANTITIES FOR THE PLOTS
 
-if nargout == 2
-    %parout.integration.time = theta;
-    
+if nargout == 2 || ~isempty(wrapper)
     parout.interp.mach = mach;
     parout.interp.alpha = alphaOut;
     parout.interp.beta = betaOut;
     parout.interp.alt = altitude;
     parout.interp.mass = m;
-    parout.interp.inertias = [Ixx, Iyy, Izz]; 
+    parout.interp.inertias = [Ixx; Iyy; Izz]; 
 
-    parout.wind.NED = [uw, vw, ww];
+    parout.wind.NED = [uw; vw; ww];
     parout.wind.body = windVels;
     
     parout.rotations.dcm = dcm;
     
     parout.velocities = vels;
     
-    parout.forces.aero = [fX, fY, fZ];
+    parout.forces.aero = [fX; fY; fZ];
     parout.forces.T = T;
     
     parout.air.rho = rho;
     parout.air.P = P;
     
-    parout.accelerations.body = [du, dv, dw];
-    parout.accelerations.angular = [dp, dq, dr];
+    parout.accelerations.body = [du; dv; dw];
+    parout.accelerations.angular = [dp; dq; dr];
     
     parout.coeff.CA = CA;
     parout.coeff.CYB = CYB;
@@ -359,5 +358,7 @@ if nargout == 2
     parout.coeff.Cnp = Cnp;
     parout.coeff.XCPlon = XCPlon;
     parout.coeff.XCPlat = XCPlat;
-    parout.coeff.XCPtot = XCPtot; 
+    parout.coeff.XCPtot = XCPtot;
+    
+    if ~isempty(wrapper), wrapper.setCache(parout); end
 end
\ No newline at end of file
diff --git a/functions/odeOutputFcn/storeData.m b/functions/odeOutputFcn/storeData.m
new file mode 100644
index 0000000000000000000000000000000000000000..54d45fb6c1b2619d2a1d103fb3462ae6f1a13c68
--- /dev/null
+++ b/functions/odeOutputFcn/storeData.m
@@ -0,0 +1,10 @@
+function s = storeData(~, ~, flag, varargin)
+s = 0;
+switch flag
+    case {'init', 'done'}
+    case ''
+        wrapper = varargin{end};
+        wrapper.setData();
+end
+end
+
diff --git a/settings/odeConfig.m b/settings/odeConfig.m
index 3d957cc27113687be806fc901f56bcdd4c82ed91..15f92ea389448bd17d66e48edb5d95950487cdbf 100644
--- a/settings/odeConfig.m
+++ b/settings/odeConfig.m
@@ -2,11 +2,13 @@
 
 % ODE settings
 ode.finalTime =  2000;                                                      % [s] Final integration time
-ode.optAscent = odeset('Events', @eventApogee, 'InitialStep', 1);           % ODE options for ascend
+ode.optAscent = odeset('Events', @eventApogee, 'InitialStep', 1, ...
+    'OutputFcn', @storeData);                                               % ODE options for ascend
 ode.optAscentDelayPara= odeset('InitialStep', 1);                           % ODE options for due to the opening delay of the parachute
-ode.optParachute = odeset('Events', @eventParaCut);                         % ODE options for the parachutes
+ode.optParachute = odeset('Events', @eventParaCut, 'InitialStep',  1, ....
+    'OutputFcn', @storeData);                                               % ODE options for the parachutes
 ode.optDescent = odeset('Events', @eventLanding,...
-'RelTol', 1e-3, 'AbsTol', 1e-3);                                            % ODE options for ballistic descent
+'RelTol', 1e-3, 'AbsTol', 1e-3, 'OutputFcn', @storeData);                   % ODE options for ballistic descent
 ode.optLaunchpad = odeset('Events', @eventPad);                             % ODE options for the launchpad phase
 
 % Settings for descent 6dof simulation