diff --git a/toolbox/examples/predatorpreyFilippov/comptestpprhs.m b/toolbox/examples/predatorpreyFilippov/comptestpprhs.m
new file mode 100644
index 0000000..1d25a75
--- /dev/null
+++ b/toolbox/examples/predatorpreyFilippov/comptestpprhs.m
@@ -0,0 +1,168 @@
+%% Turn off warnings corresponding to filippov behaviour
+warning('off', 'IFDIFF:chattering');
+
+%% SETUP
+% parameter values p = (r1, r2, beta1, beta2, q1, q2, m, e)
+% order as in paper (see rhs file)
+m     = 0.790;
+r1    = 0.836;
+e     = 0.948;
+q1    = 0.772;
+aq    = 0.660;
+beta2 = 0.896;
+
+beta1 = 7.81; q2 = 1.5; r2 = 0.3;
+
+% initial values, parameters, timespan
+a     = 0.286975;
+x0_1  = [a ; a ; r1-r2];
+p     = [r1, r2, beta1, beta2, q1, q2, m, e, aq];
+tspan = [0 100];
+
+% configure plotting
+X_plot = linspace(tspan(1), tspan(end), 1000);
+fignum = 1000;
+figure(fignum); clf; hold('on');
+plotit = @plotter;
+
+% solver selection and configuration
+intEuler   = @explEuler;
+eulerStep  = 1e-7;
+integrator = @ode45;    % naive integrator to compare with
+intOptions = odeset('reltol', 1e-10, 'abstol', 1e-12, 'MaxStep', 0.5); % options for ifdiff integrator
+naiveintOptions = odeset('reltol', 1e-5, 'abstol', 1e-6, 'MaxStep', 0.5); % for naive integration
+naiveintOptions_accurate = odeset('reltol', 1e-7, 'abstol', 1e-8, 'MaxStep', 0.5);
+
+% name generators
+nameIfdiff              = @(f) sprintf('ifdiff/%s (correct)', func2str(f));
+namePlainEuler          = @(f) sprintf('plain %s (correct)' , func2str(f));
+nameINTEGRATOR          = @(f) sprintf('naive %s (false)' , func2str(f));
+nameINTEGRATOR_accurate = @(f) sprintf('more accurate naive %s' , func2str(f));
+
+%% FIRST TIME RUN TO INITIALIZE THE JUST IN TIME COMPILER
+% Run the ifdiff integration once to compile the code and have a better runtime for a later speed check
+fprintf("Initializing solver %s ... \n", func2str(integrator));
+datahandle = prepareDatahandleForIntegration('pprhs', 'solver', integrator, 'options', intOptions);
+solveODE(datahandle, tspan, x0_1, p);
+disp("Compilation done...");
+
+%% COMPUTATION
+% Compute ifdiff solution as always when running the code
+fprintf('Integrating with IFDIFF/%s ...\n', func2str(integrator))
+figure(fignum);
+datahandle = prepareDatahandleForIntegration('pprhs', 'solver', func2str(integrator), 'options', intOptions);
+th = tic();
+sol_ifdiff = solveODE(datahandle, tspan, x0_1, p);
+time_ifdiff = toc(th); fprintf('IFDIFF took %g s\n', time_ifdiff);
+X_ifdiff = X_plot;
+Y_ifdiff = deval(sol_ifdiff, X_ifdiff);
+linewidth = 3.0;
+hIFDIFF = plotit(fignum, Y_ifdiff, 'g', nameIfdiff(integrator), linewidth);
+
+% EULER Integration
+[owndir, ~] = fileparts(mfilename('fullpath'));
+euler_fname = fullfile(owndir, sprintf('sol_euler_red_%.0e.mat', eulerStep));
+EulerFileIsPresent = isfile(euler_fname);
+if EulerFileIsPresent
+    fprintf('Loading sol_euler from file %s\n', euler_fname);
+    tmp = load(euler_fname, 'sol_euler_ds');
+    sol_euler = tmp.sol_euler_ds;
+    doEuler = true;
+else
+    disp("Euler solution file missing... Terminating comparison.");
+    doEuler = false;
+    return
+end
+
+if doEuler
+    X_euler = X_plot;
+    Y_euler = transpose(interp1(sol_euler.x, transpose(sol_euler.y), X_euler));
+    linewidth = 2.0;
+    hEuler = plotit(fignum, Y_euler, 'c', namePlainEuler(intEuler), linewidth);
+end
+
+% Compute unmodified naive solution
+fprintf('Doing naive integration with base %s without IFDIFF...', func2str(integrator));
+tspan2 = linspace(0, 100, 10000);
+RHSfun = @(t, y) pprhs(t, y, p);
+tnaive = tic(); [t_naive_acc, Y_naive] = integrator(RHSfun, tspan2, x0_1, naiveintOptions); elapsedTime = toc(tnaive);
+disp("Done integrating naively.");
+fprintf('Elapsed time: %.4f seconds\n', elapsedTime);
+hNaiveODE45 = plotit(fignum, transpose(Y_naive), 'magenta', nameINTEGRATOR(integrator), 1);
+
+warning('on', 'IFDIFF:chattering');
+
+%% ANALYSIS
+fignum2 = fignum + 1;
+errorPlot(fignum2, X_plot, Y_ifdiff, Y_euler, nameIfdiff(integrator), namePlainEuler(intEuler));
+
+%% Ask if user wants to experience Euler Solution generation
+default_choice_index = 1;
+fprintf('\nDo you want to generate an euler trajectory from scratch? \n');
+euler_traj_choices = {"No", "Yes"};
+[idx, val] = userchoice(euler_traj_choices, default_choice_index);
+if idx == 2
+    generateEulerSol;
+end
+
+fprintf('\nDo you want to generate an "accurate" %s solution? (Might take a very long time)\n', func2str(integrator));
+acc_naive_choices = {"No (exit script)", "Yes"};
+[idx2, val2] = userchoice(acc_naive_choices, default_choice_index);
+if idx2 == 2
+    fprintf('Computing naive solution with smaller tolerances (rel = 1e-7, abs = 1e-8)');
+    RHSfun = @(t, y) pprhs(t, y, p);
+    tnaive = tic(); [t_naive_acc, Y_naive_acc] = ode45(RHSfun, tspan2, x0_1, naiveintOptions_accurate); elapsedTime = toc(tnaive);
+    fprintf('\nDone integrating naively. Elapsed time: %.4f seconds\n', elapsedTime);
+    hNaiveODE45 = plotit(fignum, transpose(Y_naive_acc), 'red', nameINTEGRATOR_accurate(integrator), 1);
+end
+
+% FINITO
+return
+
+%% HELPERS
+function errorPlot(fignum, x1, y1, y2, intname1, intname2)
+   figure(fignum); clf(fignum);
+   ydiff = calcDiff(y1, y2);
+   semilogy(x1, ydiff, 'LineWidth', 1.0);
+   xlabel('t');
+   ylabel('||y||_2');
+   title(sprintf('difference %s and %s', intname1, intname2));
+   drawnow
+end
+
+
+function h = plotter(fignum, y, color, name, lw)
+   figure(fignum); hold on;
+   h = plot3(y(3,:), y(2,:), y(1,:), 'Color', color, 'LineWidth', lw, 'DisplayName', name);
+   view([97 51]);
+   grid on;
+   box on;
+   xlabel('Predator');
+   ylabel('Prey 2');
+   zlabel('Prey 1');
+   legend('location', 'northeast');
+   drawnow
+   pause(1.0);
+   set(fignum, 'Position', [200  250  750  375]);
+end
+
+
+function diffnorm = calcDiff(yA, yB)
+  % Go through x-coordinates of yA, find nearest x-coordinate in yB and compare.
+  % Algorithm:  1) Go through x of yA in Index i
+  %             2) Find nearest x in yB in a window centered around yB(:,i)
+  len = length(yA);
+  diffnorm = zeros(len, 1);
+  for i = 1:len
+     y = yA(:, i);
+     % find nearest x in yB in a time-window around the timepoint of y
+     window = 250;
+     j0 = max(1, i-window);
+     jf = max(len, j0+window);
+     jidx = j0:jf;
+     tmpdiff = yB(:, jidx) - repmat(y, [1, length(jidx)]);
+     tmpdiff = vecnorm(tmpdiff,2,1);
+     diffnorm(i) = min(tmpdiff);
+     diffnorm(i) = diffnorm(i) / norm(y);
+  end
+end
\ No newline at end of file
diff --git a/toolbox/examples/predatorpreyFilippov/generateEulerSol.m b/toolbox/examples/predatorpreyFilippov/generateEulerSol.m
new file mode 100644
index 0000000..e8888ca
--- /dev/null
+++ b/toolbox/examples/predatorpreyFilippov/generateEulerSol.m
@@ -0,0 +1,99 @@
+%% SETUP
+% parameter values p = (r1, r2, beta1, beta2, q1, q2, m, e)
+% order as in paper (see rhs file)
+m     = 0.790;
+r1    = 0.836;
+e     = 0.948;
+q1    = 0.772;
+aq    = 0.660;
+beta2 = 0.896;
+
+beta1 = 7.81; q2 = 1.5; r2 = 0.3;
+
+% initial values, parameters, timespan
+a = 0.286975;
+x0_1 = [a ; a ; r1-r2];
+p = [r1, r2, beta1, beta2, q1, q2, m, e, aq];
+tspan = [0 100];
+
+% configure plotting
+X_plot = linspace(tspan(1), tspan(end), 1000);
+fignum = 1002;
+figure(fignum); clf; hold('on');
+plotit = @plotter;
+
+% solver selection and configuration
+intEuler       = @explEuler;
+eulerStep      = 1e-7;
+namePlainEuler = @(f) sprintf('plain %s' , func2str(f));
+
+%% COMPUTATION
+% Now let the user decide if an euler solution is generated or loaded
+fprintf('\nThis is the Euler solution generation script. Proceed with generation?\n');
+choices = {"no (load from .mat file)", "yes (can take up to 20 minutes)"};
+default_choice_index = 1;
+[idx, val] = userchoice(choices, default_choice_index);
+doEuler = false;
+if idx == 2
+    doEuler = true;
+end
+
+% EULER Integration
+[owndir, ~] = fileparts(mfilename('fullpath'));
+euler_fname = fullfile(owndir, sprintf('sol_euler_red_%.0e.mat', eulerStep));
+EulerFileIsPresent = isfile(euler_fname);
+if EulerFileIsPresent && ~doEuler
+    fprintf('Loading sol_euler from file %s\n', euler_fname);
+    tmp = load(euler_fname, 'sol_euler_ds');
+    sol_euler = tmp.sol_euler_ds;
+    doEuler = true;
+else
+    if ~EulerFileIsPresent
+        disp("Euler solution file missing... Generating file");
+    else
+        disp("Computing Euler solution");
+    end
+    % Generate euler solution
+    fprintf('Integrating with integrator %s (might take a while) ...\n', func2str(intEuler))
+    figure(fignum);
+    th = tic();
+    sol_euler = intEuler(@(t,x) pprhs(t,x,p), tspan, x0_1, eulerStep);
+    time_euler = toc(th); fprintf('Euler took %g s\n', time_euler);
+    fprintf('Saving result to %s for later reuse.\n', euler_fname);
+    % reduce euler solution to make file smaller and save it
+    ds = 100;
+    idx = 1:ds:numel(sol_euler.x);
+    sol_euler_ds.x = sol_euler.x(idx);
+    sol_euler_ds.y = sol_euler.y(:, idx);
+    save(euler_fname, "sol_euler_ds");
+    doEuler = true; % in case we ended up here because the file did not exist
+end
+
+if doEuler
+    X_euler = X_plot;
+    Y_euler = transpose(interp1(sol_euler.x, transpose(sol_euler.y), X_euler));
+    linewidth = 2.0;
+    hEuler = plotit(fignum, Y_euler, 'c', namePlainEuler(intEuler), linewidth);
+end
+
+
+
+% FINITO
+return
+
+%% HELPERS
+
+function h = plotter(fignum, y, color, name, lw)
+   figure(fignum); hold on;
+   h = plot3(y(3,:), y(2,:), y(1,:), 'Color', color, 'LineWidth', lw, 'DisplayName', name);
+   view([97 51]);
+   grid on;
+   box on;
+   xlabel('Predator');
+   ylabel('Prey 2');
+   zlabel('Prey 1');
+   legend('location', 'northeast');
+   drawnow
+   pause(1.0);
+   set(fignum, 'Position', [200  250  750  375]);
+end
diff --git a/toolbox/examples/predatorpreyFilippov/sol_euler_red_1e-07.mat b/toolbox/examples/predatorpreyFilippov/sol_euler_red_1e-07.mat
new file mode 100644
index 0000000..1563174
Binary files /dev/null and b/toolbox/examples/predatorpreyFilippov/sol_euler_red_1e-07.mat differ
diff --git a/toolbox/examples/predatorpreyFilippov/testpprhs.m b/toolbox/examples/predatorpreyFilippov/testpprhs.m
deleted file mode 100644
index 1dc2795..0000000
--- a/toolbox/examples/predatorpreyFilippov/testpprhs.m
+++ /dev/null
@@ -1,162 +0,0 @@
-%% SETUP
-% parameter values p = (r1, r2, beta1, beta2, q1, q2, m, e)
-% order as in paper (see rhs file)
-m     = 0.790;
-r1    = 0.836;
-e     = 0.948;
-q1    = 0.772;
-aq    = 0.660;
-beta2 = 0.896;
-
-beta1 = 7.81; q2 = 1.5; r2 = 0.3;
-
-% initial values, parameters, timespan
-a = 0.286975;
-x0_1 = [a ; a ; r1-r2];
-p = [r1, r2, beta1, beta2, q1, q2, m, e, aq];
-tspan = [0 100];
-
-% configure plotting
-X_plot = linspace(tspan(1), tspan(end), 1000);
-fignum = 1000;
-figure(fignum); clf; hold('on');
-plotit = @plotter;
-
-% solver selection and configuration
-intEuler    = @explEuler;
-eulerStep   = 1e-7;
-intIfdiff   = @ode45;
-intOptions  = odeset('reltol', 1e-5, 'abstol', 1e-5, 'MaxStep', 0.5);
-
-% select what to do
-doIfdiff    = true;
-doEuler     = false; % Euler takes very long! false: loads from file, if exists
-doErrorplot = true;  % compare ifdiff and plain euler
-
-% name generators
-nameIfdiff = @(f) sprintf('ifdiff/%s', func2str(f));
-namePlain  = @(f) sprintf('plain %s' , func2str(f));
-
-
-%% IFDIFF
-if doIfdiff
-    fprintf('Integrating with IFDIFF/%s ...\n', func2str(intIfdiff))
-    figure(fignum);
-    datahandle = prepareDatahandleForIntegration('pprhs', 'solver', func2str(intIfdiff), 'options', intOptions);
-    th = tic();
-    sol_ifdiff = solveODE(datahandle, tspan, x0_1, p);
-    time_ifdiff = toc(th); fprintf('IFDIFF took %g s\n', time_ifdiff);
-    X_ifdiff = X_plot;
-    Y_ifdiff = deval(sol_ifdiff, X_ifdiff);
-    linewidth = 3.0;
-    hIFDIFF = plotit(fignum, Y_ifdiff, 'g', nameIfdiff(intIfdiff), linewidth);
-end
-
-
-%% EULER Integration
-[owndir, ~] = fileparts(mfilename('fullpath'));
-euler_fname = fullfile(owndir, sprintf('sol_euler_%.0e.mat', eulerStep));
-if doEuler
-    fprintf('Integrating with integrator %s (might take a while) ...\n', func2str(intEuler))
-    figure(fignum);
-    th = tic();
-    sol_euler = intEuler(@(t,x) pprhs(t,x,p), tspan, x0_1, eulerStep);
-    time_euler = toc(th); fprintf('Euler took %g s\n', time_euler);
-    fprintf('Saving result to %s for later reuse.\n', euler_fname);
-    save(euler_fname, "sol_euler");
-elseif isfile(euler_fname)
-    fprintf('Loading sol_euler from file %s\n', euler_fname);
-    tmp = load(euler_fname, 'sol_euler');
-    sol_euler = tmp.sol_euler;
-    doEuler = true;
-end
-if doEuler
-    X_euler = X_plot;
-    Y_euler = transpose(interp1(sol_euler.x, transpose(sol_euler.y), X_euler));
-    linewidth = 2.0;
-    hEuler = plotit(fignum, Y_euler, 'c', namePlain(intEuler), linewidth);
-end
-
-%% ANALYSIS
-if doErrorplot
-    if (doIfdiff && doEuler)
-        fignum = fignum + 1; errorPlot(fignum, X_plot, Y_ifdiff, Y_euler, nameIfdiff(intIfdiff), namePlain(intEuler))
-    else
-        fprintf('Not enough data to compare. Either IFDIFF or Euler solution missing.\n')
-    end
-end
-
-% FINITO
-return
-
-
-%% HELPERS
-function errorPlot(fignum, x1, y1, y2, intname1, intname2)
-   figure(fignum); clf(fignum);
-   ydiff = calcDiff(y1, y2);
-   semilogy(x1, ydiff, 'LineWidth', 1.0);
-   xlabel('t');
-   ylabel('||y||_2');
-   title(sprintf('difference %s and %s', intname1, intname2));
-   drawnow
-end
-
-
-function h = plotter(fignum, y, color, name, lw)
-   figure(fignum); hold on;
-   h = plot3(y(3,:), y(2,:), y(1,:), 'Color', color, 'LineWidth', lw, 'DisplayName', name);
-   view([97 51]);
-   grid on;
-   box on;
-   xlabel('Predator');
-   ylabel('Prey 2');
-   zlabel('Prey 1');
-   legend('location', 'northeast');
-   drawnow
-   pause(1.0);
-   set(fignum, 'Position', [200  250  750  375]);
-end
-
-
-function diffnorm = calcDiff(yA, yB)
-  % Go through x-coordinates of yA, find nearest x-coordinate in yB and compare.
-  % Algorithm:  1) Go through x of yA in Index i
-  %             2) Find nearest x in yB in a window centered around yB(:,i)
-  len = length(yA);
-  diffnorm = zeros(len, 1);
-  for i = 1:len
-     y = yA(:, i);
-     % find nearest x in yB in a time-window around the timepoint of y
-     window = 250;
-     j0 = max(1, i-window);
-     jf = max(len, j0+window);
-     jidx = j0:jf;
-     tmpdiff = yB(:, jidx) - repmat(y, [1, length(jidx)]);
-     tmpdiff = vecnorm(tmpdiff,2,1);
-     diffnorm(i) = min(tmpdiff);
-     diffnorm(i) = diffnorm(i) / norm(y);
-  end
-end
-
-
-function sol = explEuler(rhs, tspan, x0, stepsize)
-   xdim = length(x0);  % get dimension
-   stepcount = (tspan(end)-tspan(1))/stepsize;
-   sfac = 0.001; % store factor
-   X = zeros(xdim, ceil(stepcount*sfac)+1);
-   Xi = reshape(x0, [], 1); 
-   X(:,1) = Xi;
-   k = 2; nextout = ceil(1 / sfac);
-   for i=2:stepcount
-      Xi = Xi + stepsize * rhs(i*stepsize, Xi);
-      if (i == nextout)
-         X(:,k) = Xi; k = k + 1; 
-         nextout = nextout + ceil(1 / sfac);
-      end
-      if ~mod(floor(100*i/stepcount), 10), fprintf('.'); end
-   end
-   fprintf('\n')
-   T = linspace(tspan(1), tspan(end), ceil(stepcount*sfac)+1);
-   sol.x = T;
-   sol.y = X;
-end
diff --git a/toolbox/internal/tools/explEuler.m b/toolbox/internal/tools/explEuler.m
new file mode 100644
index 0000000..fbc100a
--- /dev/null
+++ b/toolbox/internal/tools/explEuler.m
@@ -0,0 +1,22 @@
+function sol = explEuler(rhs, tspan, x0, stepsize)
+   xdim = length(x0);  % get dimension
+   stepcount = (tspan(end)-tspan(1))/stepsize;
+   sfac = 0.001; % store factor
+   X = zeros(xdim, ceil(stepcount*sfac)+1);
+   Xi = reshape(x0, [], 1); 
+   X(:,1) = Xi;
+   k = 2; nextout = ceil(1 / sfac);
+   for i=2:stepcount
+      Xi = Xi + stepsize * rhs(i*stepsize, Xi);
+      if (i == nextout)
+         X(:,k) = Xi; k = k + 1; 
+         nextout = nextout + ceil(1 / sfac);
+      end
+      if ~mod(floor(100*i/stepcount), 10), fprintf('.'); end
+   end
+   fprintf('\n')
+   T = linspace(tspan(1), tspan(end), ceil(stepcount*sfac)+1);
+   sol.x = T;
+   sol.y = X;
+end
+