diff --git a/.gitignore b/.gitignore
index e6907a0..d6c8d60 100644
--- a/.gitignore
+++ b/.gitignore
@@ -35,6 +35,7 @@
 !/test/fnft__poly/
 !/test/fnft_version_test/
 !/test/fnft_version_test.c/
+!/test/fnft_kdvv_inverse/
 
 # Ignore MacOS hidden files in subdirectories
 /**/*.DS_Store
diff --git a/examples/fnft_kdvv_inverse_example_1.c b/examples/fnft_kdvv_inverse_example_1.c
new file mode 100644
index 0000000..52e1d26
--- /dev/null
+++ b/examples/fnft_kdvv_inverse_example_1.c
@@ -0,0 +1,75 @@
+/*
+* This file is part of FNFT.  
+*                                                                  
+* FNFT is free software; you can redistribute it and/or
+* modify it under the terms of the version 2 of the GNU General
+* Public License as published by the Free Software Foundation.
+*
+* FNFT is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*                                                                      
+* You should have received a copy of the GNU General Public License
+* along with this program. If not, see <http://www.gnu.org/licenses/>.
+*
+* Contributors:
+* Fabian Fischer (Hiwi KIT) 2026
+*/
+
+// This is a simple example for the usage of the inverse kdvv transform in C
+
+#define FNFT_ENABLE_SHORT_NAMES
+
+#include <stdio.h>
+
+#include "fnft_kdvv_inverse.h"
+#include "fnft__errwarn.h"
+#include "fnft__misc.h"
+
+#define K 5     // number of bound states to add
+
+INT main()
+{
+    INT ret_code = SUCCESS;
+    
+    COMPLEX * q = NULL;
+    COMPLEX * contspec = NULL;
+
+    // five desired bound states and norming constants
+    // - bound states have to be positive, purely imaginary numbers
+    // - bound states have to be in descending order
+    // - every norming constants belong to the bound state at the same index
+    // - the signs of the norming constants have to alternate regarding the order
+    //   of bound states. The norming constant of the biggest bound state has to be positive.
+    COMPLEX bound_states[K] = { I*5.0, I*4.0, I*3.0, I*2.0, I*1.0 };
+    COMPLEX norming_constants[K] = { 1e5, -1e-2, 1e0, -1e3, 1e1 };
+
+    // General parameters
+    UINT D = 256;                   // Number of samples of the computed output
+    REAL T[2] = {-10.0, 10.0};      // area for which the output should be computed
+    
+    // allocation of memory for the computed output
+    q = malloc(D * sizeof(COMPLEX));
+    CHECK_NOMEM(q, ret_code, leave_fun);
+    
+    // Define continuous spectrum
+    // continuous spectrum is out of function, but needs to be defined/allocated for using 
+    // the inverse kdvv (state 04/2026)
+    UINT M = 10;
+    REAL XI[2] = {-2.0, 2.0};
+
+    contspec = malloc(M * sizeof(COMPLEX));
+    CHECK_NOMEM(contspec, ret_code, leave_fun);
+
+    // call of the inverse kdvv
+    ret_code = fnft_kdvv_inverse(M, contspec, XI, K, bound_states, norming_constants, D, q, T, NULL);
+    CHECK_RETCODE(ret_code, leave_fun);
+
+    // prints the results in the console
+    misc_print_buf(D, q, "output");
+
+leave_fun:
+    free(q);
+    free(contspec);
+}
diff --git a/examples/fnft_kdvv_inverse_example_2.c b/examples/fnft_kdvv_inverse_example_2.c
new file mode 100644
index 0000000..04de867
--- /dev/null
+++ b/examples/fnft_kdvv_inverse_example_2.c
@@ -0,0 +1,88 @@
+/*
+* This file is part of FNFT.  
+*                                                                  
+* FNFT is free software; you can redistribute it and/or
+* modify it under the terms of the version 2 of the GNU General
+* Public License as published by the Free Software Foundation.
+*
+* FNFT is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*                                                                      
+* You should have received a copy of the GNU General Public License
+* along with this program. If not, see <http://www.gnu.org/licenses/>.
+*
+* Contributors:
+* Fabian Fischer (Hiwi KIT) 2026
+*/
+
+// This is an example for the usage of the inverse kdvv transform with many 
+// bound states
+
+#define FNFT_ENABLE_SHORT_NAMES
+
+#include <stdio.h>
+
+#include "fnft_kdvv_inverse.h"
+#include "fnft__errwarn.h"
+#include "fnft__misc.h"
+
+#define K 19    // Number of bound states to add
+
+INT main()
+{
+    INT ret_code = SUCCESS;
+
+    COMPLEX * bound_states = NULL;
+    COMPLEX * q = NULL;
+    COMPLEX * contspec = NULL;
+
+    COMPLEX desired_solitions_height[K] = { 40.0, 39.0, 38.0, 37.0, 36.0,
+                                            30.0, 29.0, 28.0, 27.0, 26.0,
+                                            20.0, 19.0, 18.0, 17.0, 16.0,
+                                            10.0, 9.0, 8.0, 7.0};
+
+    COMPLEX normconsts[K] = {   1e20, -1e-7, 1e5, -1e3, 1e1,
+                                -1e0, 1e2, -1e4, 1e-6, -1e8,
+                                1e2, -1e4, 1e6, -1e8, 1e-10,
+                                -1e7, 1e-6, -1e5, 1e-9};
+
+    // resulting bound states out of desired solitions height                                
+    bound_states = malloc(K * sizeof(COMPLEX));
+    CHECK_NOMEM(bound_states, ret_code, leave_fun);
+                                
+    for (UINT i = 0; i<K; i++){
+        bound_states[i] = I*SQRT(desired_solitions_height[i]/2.0);
+    }
+
+    // General parameters
+    UINT D = 256;
+    REAL T[2] = {-18.0, 12.0};
+    
+    // allocating memory for the output of the inverse kdvv
+    q = malloc(D * sizeof(COMPLEX));
+    CHECK_NOMEM(q, ret_code, leave_fun);
+    
+    // Define Parameters and allocate memory for continuous spectrum
+    // inverse kdvv cannot deal with a nonzero continuous spectrum (state 04/2026), 
+    // hence it is implicitly assumed to zero.
+    // Nevertheless the variables needs to be delivered
+    UINT M = 10;
+    REAL XI[2] = {-2.0, 2.0};
+
+    contspec = malloc(M * sizeof(COMPLEX));
+    CHECK_NOMEM(contspec, ret_code, leave_fun);
+
+    // call inverse kdvv transform
+    ret_code = fnft_kdvv_inverse(M, contspec, XI, K, bound_states, normconsts, D, q, T, NULL);
+    CHECK_RETCODE(ret_code, leave_fun);
+
+    // print output of inverse kdvv into console
+    misc_print_buf(D, q, "output");
+
+leave_fun:
+    free(bound_states);
+    free(q);
+    free(contspec);
+}
diff --git a/examples/mex_fnft_kdvv_inverse_example_1.m b/examples/mex_fnft_kdvv_inverse_example_1.m
new file mode 100644
index 0000000..6bb45f9
--- /dev/null
+++ b/examples/mex_fnft_kdvv_inverse_example_1.m
@@ -0,0 +1,73 @@
+% This file is part of FNFT.
+%
+% FNFT is free software; you can redistribute it and/or
+% modify it under the terms of the version 2 of the GNU General
+% Public License as published by the Free Software Foundation.
+%
+% FNFT is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+% Contributors:
+% Sander Wahls (KIT) 2026.
+% Fabian Fischer (Hiwi KIT) 2026.
+
+% This example should illustrate how to use the inverse kdvv transform
+
+clear all;
+close all;
+
+% desired height of solitions:
+desired_solitions = [4, 3, 2, 1];
+
+% resulting bound states out of desired heights of solitions:
+% - bound states have to be positive, purely imaginary numbers
+% - the bound states have to be in descendent order
+bound_states = 1i*sqrt(desired_solitions ./2);
+
+% desired norming constants
+% defines how much the solitions are shifted towards each other
+% - signs have to alternate regards to the order of the bound states
+% - sign of the normconst for the biggest eigenvalue has to be positive
+norming_constants = complex([1, -1, 1, -1].*[0.00001, 0.1, 1, 10]);
+
+% Number of samples of the output of the inverse kdvv
+D = 1001;
+
+% Area, for which the output has to be computed
+% Can be asymmetric
+T = [-10 10];
+
+% defining the continuous spectrum:
+% out of function, just for seek of completeness (state 04/2026)
+contspec = [];
+XI = [1e-6 1];
+
+% calls function of c-library FNFT
+q = mex_fnft_kdvv_inverse(contspec, XI, bound_states, norming_constants, D, T);
+% the result consists of solitions with desired height (if the solitions are far
+% away enough to each other)
+
+% --- Plot the results ---
+t = linspace(T(1), T(2), D);
+
+% the output of the inverse kdvv is defined as a complex number (state 04/2026),
+% however only the imaginary part should be zero and can be neglected
+q = double(real(q(:)'));
+
+figure;
+plot(t, q);
+title('output of inverse kdvv: time-domain');
+xlabel('t');
+ylabel('q(t)');
+legend('q(t)');
+
+figure;
+stem(imag(bound_states),real(norming_constants), 'x');
+title('Bound states and norming constants');
+xlabel('bound states');
+ylabel('norming constants');
\ No newline at end of file
diff --git a/examples/mex_fnft_kdvv_inverse_example_2.m b/examples/mex_fnft_kdvv_inverse_example_2.m
new file mode 100644
index 0000000..fcc8811
--- /dev/null
+++ b/examples/mex_fnft_kdvv_inverse_example_2.m
@@ -0,0 +1,85 @@
+% This file is part of FNFT.
+%
+% FNFT is free software; you can redistribute it and/or
+% modify it under the terms of the version 2 of the GNU General
+% Public License as published by the Free Software Foundation.
+%
+% FNFT is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+% Contributors:
+% Sander Wahls (KIT) 2026.
+% Fabian Fischer (Hiwi KIT) 2026.
+
+% This example should illustrate, how to use the inverse kdvv in combination 
+% with the forward kdvv transform. For basic usage of the inverse kdvv, please
+% please have a look at mex_fnft_kdvv_inverse_example_1.m (state 04/2026).
+
+clear all;
+close all;
+
+contspec_initial = [];
+XI = [1e-6 10];
+
+bound_states_initial = 1i*sqrt([4, 3, 2, 1] ./2);
+norming_constants_initial = complex([1, -1, 1, -1].*[10000, 0.1, 1, 0.001]);
+
+D = 1001;
+T = [-15 10];
+
+% inverse kdvv transform
+q = mex_fnft_kdvv_inverse(contspec_initial, XI, bound_states_initial, norming_constants_initial, D, T);
+
+q = double(real(q(:)'));
+
+% compute the nonlinear Fourier transform of the output of the inverse kdvv
+[contspec_computed, bound_states_computed, norming_constants_computed] = mex_fnft_kdvv(q, T, XI);
+
+
+% --- Plot the results ---
+t = linspace(T(1), T(2), D);
+
+figure;
+plot(t, q);
+title('output of inverse kdvv: time-domain');
+xlabel('t');
+ylabel('q(t)');
+legend('q(t)');
+
+% plotting the initial and the computed discrete spectrum
+% the difference between the initial and the computed spectrum should be very small
+
+figure;
+stem(imag(bound_states_initial),real(norming_constants_initial), 'x', 'color', 'r');
+hold on;
+stem(imag(bound_states_computed),real(norming_constants_computed), 'x', 'color', 'g');
+hold off;
+title('initial and computed bound states and norming constants');
+xlabel('bound states');
+ylabel('norming constants');
+legend('initial', 'computed');
+
+% plotting of continuous spectrum
+% values are small because initial continuous spectrum was assumed as zero
+% and the computed continuous spectrum is only the result of numerical errors
+
+ep_xi = (XI(2) - XI(1)) / (D - 1);
+xi = XI(1):ep_xi:XI(2);
+
+figure;
+hmag=subplot(2,1,1);
+semilogy(xi, abs(contspec_computed));
+title('computed continuous spectrum');
+xlabel('\xi');
+ylabel('|r(\xi)|');
+hang=subplot(2,1,2);
+plot(xi, angle(contspec_computed));
+xlabel('\xi');
+ylabel('\angle r(\xi)');
+linkaxes([hmag,hang],'x');
+
diff --git a/examples/mex_fnft_kdvv_inverse_example_3.m b/examples/mex_fnft_kdvv_inverse_example_3.m
new file mode 100644
index 0000000..9334d3a
--- /dev/null
+++ b/examples/mex_fnft_kdvv_inverse_example_3.m
@@ -0,0 +1,78 @@
+% This file is part of FNFT.
+%
+% FNFT is free software; you can redistribute it and/or
+% modify it under the terms of the version 2 of the GNU General
+% Public License as published by the Free Software Foundation.
+%
+% FNFT is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+% Contributors:
+% Fabian Fischer (Hiwi KIT) 2026.
+
+% This example shows a fast way to try different values (bound states, norming constants)
+% for the inverse kdvv transform, because this script sorts and considers the alternating 
+% signs of the norming constants automatically.
+
+clear all;
+close all;
+
+% === modify here ====================================================================================
+% a solition in the output of inverse kdvv is defined by the elements of both arrays at the same index
+desired_solitions = [9, 2, 4, 3, 7];
+desired_shifting = [10000, 0.1, 1, 0.00001, 2];
+% ====================================================================================================
+
+
+% Check if both arrays contains the same number of elements
+if length(desired_shifting) ~= length(desired_solitions)
+    error('Number of desired solitions is not equal to number of values given for shifting!');
+end
+
+% --- determining sorted bound states and norming constants out of desired values ---
+
+number_values = length(desired_solitions);
+
+tmp_bound_states = 1i*sqrt(desired_solitions ./2);
+
+% sorting
+[bound_states, indices] = sort(tmp_bound_states, "descend");
+tmp_norming_constants = desired_shifting(indices);
+
+% alternating sequence 1 and -1 as elements and with 1 as first element
+tmp_seq = (-1).^( 0:(number_values-1) );
+
+norming_constants = complex(tmp_seq.*tmp_norming_constants);
+
+% inverse kdvv
+contspec = [];
+XI = [1e-6 1];
+D = 1001;
+T = [-12 8];
+q = mex_fnft_kdvv_inverse(contspec, XI, bound_states, norming_constants, D, T);
+
+
+% --- Plot the results ---
+t = linspace(T(1), T(2), D);
+
+% the output of the inverse kdvv is defined as a complex number (state 04/2026),
+% however only the imaginary part should be zero and can be neglected
+q = double(real(q(:)'));
+
+figure;
+plot(t, q);
+title('output of inverse kdvv: time-domain');
+xlabel('t');
+ylabel('q(t)');
+legend('q(t)');
+
+figure;
+stem(imag(bound_states),real(norming_constants), 'x');
+title('Bound states and norming constants');
+xlabel('bound states');
+ylabel('norming constants');
\ No newline at end of file
diff --git a/examples/mex_fnft_kdvv_inverse_example_4.m b/examples/mex_fnft_kdvv_inverse_example_4.m
new file mode 100644
index 0000000..0ca29a4
--- /dev/null
+++ b/examples/mex_fnft_kdvv_inverse_example_4.m
@@ -0,0 +1,37 @@
+%% 
+% This file is part of FNFT.
+%
+% FNFT is free software; you can redistribute it and/or
+% modify it under the terms of the version 2 of the GNU General
+% Public License as published by the Free Software Foundation.
+%
+% FNFT is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+% Contributors:
+% Sander Wahls (KIT) 2026.
+
+% This is a simple example of how to use the inverse kdvv.
+
+close all;
+
+contspec = [];
+XI = [0 1];
+
+bound_states = 1i*sqrt( [5, 4, 3, 2, 1] /2); 
+
+normconsts = complex([1, -1, 1, -1, 1].*[10000000, 0.1, 1, 0.001, 10]);
+
+D = 1001;
+T = [-10 10];
+q = mex_fnft_kdvv_inverse(contspec, XI, bound_states, normconsts, D, T);
+
+t = linspace(T(1), T(2), D);
+plot(t, q)
+xlabel('t')
+ylabel('q(t)')
diff --git a/include/fnft_kdvv.h b/include/fnft_kdvv.h
index 542ac97..cedb0f2 100644
--- a/include/fnft_kdvv.h
+++ b/include/fnft_kdvv.h
@@ -270,7 +270,7 @@ fnft_kdvv_opts_t fnft_kdvv_default_opts();
  *       - fnft_kdv_discretization_TES4(_VANILLA)
  *
  * The accuray of the computed quantities for a given signal depends primarily on the number of samples \f$ D\f$ and the numerical method. When the exact spectrum is
- * is know, the accuracy can be quantified by defining a suitable error. The error usually decreases with increasing \f$ D\f$ assuming everthing else remains the same.
+ * known, the accuracy can be quantified by defining a suitable error. The error usually decreases with increasing \f$ D\f$ assuming everthing else remains the same.
  * The rate at which the error decreases with increase in \f$ D\f$ is known as the order of the method. The orders of the various discretizations can be found at \link fnft_kdv_discretization_t \endlink.
  * The orders of the discretizations which use exponential splitting schemes should be the same as their base methods but can deviate when accuracy of the splitting scheme is low.
  *
diff --git a/include/fnft_kdvv_inverse.h b/include/fnft_kdvv_inverse.h
new file mode 100644
index 0000000..b566788
--- /dev/null
+++ b/include/fnft_kdvv_inverse.h
@@ -0,0 +1,106 @@
+/*
+ * This file is part of FNFT.
+ *
+ * FNFT is free software; you can redistribute it and/or
+ * modify it under the terms of the version 2 of the GNU General
+ * Public License as published by the Free Software Foundation.
+ *
+ * FNFT is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ * Contributors:
+ * Sander Wahls (KIT) 2026
+ * Fabian Fischer (Hiwi KIT) 2026
+ */
+
+/**
+ * @file fnft_kdvv_inverse.h
+ * @brief Fast inverse nonlinear Fourier transform for the vanishing
+ *  Korteweg-de Vries equation.
+ * @ingroup fnft_inverse
+ */
+
+#ifndef FNFT_KDVV_INVERSE_H
+#define FNFT_KDVV_INVERSE_H
+
+#include "fnft__errwarn.h"
+#include "fnft__misc.h"
+#include "fnft_numtypes.h"
+
+/**
+ * @brief Fast inverse nonlinear Fourier transform for the 
+ *  Korteweg-de Vries equation with vanishing boundary conditions.
+ * 
+ * This routine computes the inverse nonlinear Fourier transform for the
+ * Korteweg-de Vries equation
+ * \f[ q_x + 6qq_{t} + q_{ttt}=0, \quad  q=q(x,t), \f]
+ * of Gardner et al. (<a href="https://doi.org/10.1103/PhysRevLett.19.1095">
+ * Phys. Rev. Lett., 1967</a>)
+ * 
+ * The Fast inverse nonlinear Fourier transform for the Korteweg-de Vries equation
+ * uses the crum transformation. The main references are:
+ *      - Prins and Wahls, <a href="https://doi.org/10.1016/j.cnsns.2021.105782">&quot;An accurate O(N^2) floating point algorithm for the Crum transform of the KdV equation,&quot;</a> Communications in Nonlinear Science and Numerical Simulation 102, Article 105782, 2021.
+ *      - The matlab project of P. Prins to the Crum Transformation
+ *
+ * @param[in] M Number of samples of the continuous spectrum.
+ * @param[in,out] contspec Array of length M, contains samples
+ *  \f$ \hat{q}(\xi_n) \f$, where \f$ \xi_n = XI[0] + n(XI[1]-XI[0])/(M-1) \f$
+ *  and \f$n=0,1,\dots,M-1\f$, of the to-be-inverted continuous spectrum in
+ *  ascending order (i.e.,
+ *  \f$ \hat{q}(\xi_0), \hat{q}(\xi_1), \dots, \hat{q}(\xi_{M-1}) \f$).
+ *  Note: contspec functionality currently out of function (state 04/2026)! 
+ *  It will be implicitly assumed to \f$ 0 \f$ for all \f$ \xi \f$.
+ * @param[in] XI Array of length 2, contains the position of the first and the last
+ *  sample of the continuous spectrum.
+ *  Note: contspec functionality currently out of function (state 04/2026)! It does not
+ *  matter which values are chosen at the moment.
+ * @param[in] K Number of discrete spectrum points.
+ * @param[in] bound_states Complex array of length K. Bound states have to be positive, 
+ *  purely imaginary numbers (lie on the upper half of the imaginary axis). The bound 
+ *  states have to be in descending order. To add a solition with height \f$ h_i \f$ 
+ *  the bound state have to be \f$ \gamma_i = \sqrt{ h_i/2} \f$. 
+ * @param[in] norming_constants Complex array of length K. Values of
+ *  either the norming constants \f$ b(\xi) \f$ or the residues
+ *  \f$ \frac{b(\xi)}{\partial{a(\xi)}/\partial{\xi}}\f$ at the values bound_states.
+ *  The signs of the norming constants have to alternate regards to the order of the 
+ *  bound states. The sign of the normconst for the biggest eigenvalue has to be positive
+ *  Note: currently this array is always interpreted as norming constants (state 04/2026)! 
+ *  Residues functionality is not implemented yet!
+ * @param[in] D Number of samples of the to be generated signal q.
+ * @param[out] q Array of length D. Is filled with samples
+ *  \f$ q(t_n) \f$, where \f$ t_n = T[0] + n(T[1]-T[0])/(D-1) \f$
+ *  and \f$n=0,1,\dots,D-1\f$, of the to-be-generated signal in ascending order
+ *  (i.e., \f$ q(t_0), q(t_1), \dots, q(t_{D-1}) \f$).
+ *  Has to be preallocated by the user.
+ * @param[in] T Array of length 2, contains the position in time of the first and
+ *  of the last sample of q. It should be \f$ T[0]<T[1] \f$.
+ * @param[in] opts_ptr Pointer to a \link fnft_kdvv_inverse_opts_t \endlink
+ *  object. The object  can be used to modify the behavior of the routine. Use
+ *  the routine \link fnft_kdvv_inverse_default_opts \endlink
+ *  to generate such an object and modify as desired. It is also possible to
+ *  pass NULL, in which case the routine will use the default options. The
+ *  user is reponsible to freeing the object after the routine has returned.
+ *  Note: opts_ptr is currently out of function (state 04/2026)!
+ * @return \link FNFT_SUCCESS \endlink or one of the FNFT_EC_... error codes
+ *  defined in \link fnft_errwarn.h \endlink.
+ *
+ * @ingroup fnft_inverse
+ */
+FNFT_INT fnft_kdvv_inverse(
+    const FNFT_UINT M,
+    FNFT_COMPLEX * const contspec,
+    FNFT_REAL const * const XI,
+    FNFT_UINT const K,
+    FNFT_COMPLEX const * const bound_states,
+    FNFT_COMPLEX const * const norming_constants,
+    const FNFT_UINT D,
+    FNFT_COMPLEX * const q,
+    FNFT_REAL const * const T,
+    void *opts_ptr);
+
+#endif
diff --git a/include/private/fnft__errwarn.h b/include/private/fnft__errwarn.h
index 1bf1e80..9c5643e 100644
--- a/include/private/fnft__errwarn.h
+++ b/include/private/fnft__errwarn.h
@@ -47,6 +47,12 @@
  */
 #define FNFT__E_INVALID_ARGUMENT(name) FNFT__ERRMSG(FNFT_EC_INVALID_ARGUMENT, FNFT__E_INVALID_ARGUMENT_(name))
 
+/**
+ * Macro for raising an invalid argument error with detailed description.
+ * @ingroup private_errwarn
+ */
+#define FNFT__E_INVALID_ARGUMENT_MSG(name,msg) FNFT__ERRMSG(FNFT_EC_INVALID_ARGUMENT, FNFT__E_INVALID_ARGUMENT_MSG_(name, msg))
+
 /**
  * Macro for raising an subroutine error, where ec is the error code returned
  * by the subroutine. The error code returned by this macro is -abs(ec). Since
@@ -111,6 +117,7 @@
 #define WARN(msg)           FNFT__WARN(msg)
 #define E_NOMEM             FNFT__E_NOMEM
 #define E_INVALID_ARGUMENT(name) FNFT__E_INVALID_ARGUMENT(name)
+#define E_INVALID_ARGUMENT_MSG(name,msg) FNFT__E_INVALID_ARGUMENT_MSG(name, msg)
 #define E_SUBROUTINE(ec)    FNFT__E_SUBROUTINE(ec)
 #define E_DIV_BY_ZERO       FNFT__E_DIV_BY_ZERO
 #define E_TEST_FAILED       FNFT__E_TEST_FAILED
@@ -139,6 +146,13 @@
  */
 #define FNFT__E_INVALID_ARGUMENT_(name) "Invalid argument "#name"."
 
+/**
+ * Auxiliary macro that is used to stringify the input to
+ * \link FNFT__E_INVALID_ARGUMENT_MSG \endlink. Do not call directly.
+ * @ingroup private_errwarn
+ */
+#define FNFT__E_INVALID_ARGUMENT_MSG_(name,msg) "Invalid argument "#name". "#msg
+
 /**
  * Auxiliary macro used to stringify the input to
  * \link FNFT__E_NOT_YET_IMPLEMENTED \endlink. Do not call directly.
diff --git a/include/private/fnft__kdvv_inverse_testcases.h b/include/private/fnft__kdvv_inverse_testcases.h
new file mode 100644
index 0000000..db131c8
--- /dev/null
+++ b/include/private/fnft__kdvv_inverse_testcases.h
@@ -0,0 +1,97 @@
+/*
+* This file is part of FNFT.  
+*                                                                  
+* FNFT is free software; you can redistribute it and/or
+* modify it under the terms of the version 2 of the GNU General
+* Public License as published by the Free Software Foundation.
+*
+* FNFT is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*                                                                      
+* You should have received a copy of the GNU General Public License
+* along with this program. If not, see <http://www.gnu.org/licenses/>.
+*
+* Contributors:
+* Sander Wahls (TU Delft) 2018.
+*/
+
+#ifndef FNFT__KDVV_INVERSE_TESTCASES_H
+#define FNFT__KDVV_INVERSE_TESTCASES_H
+
+#include "fnft_kdvv_inverse.h"
+#include "fnft_kdvv.h"
+
+/**
+ * @struct fnft_kdvv_params
+ */
+typedef struct {
+    UINT D;
+    REAL T[2];
+    UINT K;
+    COMPLEX * bound_states;
+    COMPLEX * normconsts;
+    UINT M;
+    REAL XI[2];
+    COMPLEX * contspec;
+} fnft_kdvv_params;
+
+
+/**
+ * @brief Routine to run tests for \link fnft_kdvv_inverse \endlink.
+ *
+ * This routine is used by the tests for \link fnft_kdvv_inverse \endlink.
+ *
+ * @param[in] params_i \link fnft_kdvv_params \endlink
+ * @param[in] err_bnd_bound_states 
+ * @param[in] err_bnd_spurious_bound_states 
+ * @param[in] err_bnd_normconst
+ * @param[in] err_bnd_contspec
+ * @return If all errors stay below bounds the routine
+ * \link FNFT_SUCCESS \endlink. Otherwise, it returns an error code
+ * (normally, \link FNFT_EC_TEST_FAILED \endlink).
+ *
+ * @ingroup kdv
+ */
+FNFT_INT fnft__kdvv_inverse_testcases_get_spectrum_of_inverse(
+    const fnft_kdvv_params params_i,
+    const FNFT_REAL err_bnd_bound_states,
+    const FNFT_REAL err_bnd_spurious_bound_states,
+    const FNFT_REAL err_bnd_normconst,
+    const FNFT_REAL err_bnd_contspec);
+ 
+    
+/**
+ * @brief Routine to print the continuous spectrum and discrete spectrum, consisting of 
+ * bound states and norming constants, as a result of \link fnft_kdvv \endlink.
+ *
+ * This routine is used by the tests for \link fnft_kdvv_inverse \endlink.
+ *
+ * @param[in] bound_states
+ * @param[in] normconsts 
+ * @param[in] contspec
+ * @param[in] XI Array of length 2, contains the position of the first and the last
+ *  sample of the continuous spectrum.
+ * @param[in] M Number of points at which the continuous spectrum is computed.
+ * @param[in] D Number of samples of the potential.
+ * @param[in] K Number of bound states (same than number of norming constants)
+ * @return void
+ *
+ * @ingroup kdv
+ */      
+void fnft__kdvv_print_spectrum(  FNFT_COMPLEX const * const bound_states,
+                            FNFT_COMPLEX const * const normconsts,
+                            FNFT_COMPLEX const * const contspec,
+                            FNFT_REAL const * const XI,
+                            const UINT M,
+                            const UINT D,
+                            const UINT K);
+
+
+#ifdef FNFT_ENABLE_SHORT_NAMES
+#define kdvv_testcases_get_spectrum_of_inverse(...) fnft__kdvv_inverse_testcases_get_spectrum_of_inverse(__VA_ARGS__)
+#define kdvv_print_spectrum(...) fnft__kdvv_print_spectrum(__VA_ARGS__)
+#endif
+
+#endif
diff --git a/include/private/fnft__misc.h b/include/private/fnft__misc.h
index 5dc9b2a..7a3ae9f 100644
--- a/include/private/fnft__misc.h
+++ b/include/private/fnft__misc.h
@@ -98,6 +98,25 @@ FNFT_REAL fnft__misc_hausdorff_dist(const FNFT_UINT lenA,
     FNFT_COMPLEX const * const vecA, const FNFT_UINT lenB,
     FNFT_COMPLEX const * const vecB);
 
+
+/**
+ * @brief Hausdorff distance between two vectors, normed by the considered element.
+ *
+ * @ingroup misc
+ * This function computes the Hausdorff distance between two vectors vecA and vecB. 
+ * By calculating the maximum distance, the distance is divided by the absolute of the 
+ * according element.
+ * @param[in] lenA Length of vector vecA.
+ * @param[in] vecA Complex vector of length lenA.
+ * @param[in] lenB length of vector vecB.
+ * @param[in] vecB Complex vector of length lenB.
+ * @return Returns the real valued Hausdorff distance between the vectors vecA and vecB.
+ */
+FNFT_REAL fnft__misc_hausdorff_dist_normed(const FNFT_UINT lenA,
+    FNFT_COMPLEX const * const vecA, const FNFT_UINT lenB,
+    FNFT_COMPLEX const * const vecB);
+
+
 /**
  * @brief Hyperbolic secant.
  *
@@ -453,6 +472,7 @@ static inline FNFT_INT fnft__misc_normalize_vector(const FNFT_UINT len, FNFT_COM
 #define misc_rel_err(...) fnft__misc_rel_err(__VA_ARGS__)
 #define misc_rel_err_real(...) fnft__misc_rel_err_real(__VA_ARGS__)
 #define misc_hausdorff_dist(...) fnft__misc_hausdorff_dist(__VA_ARGS__)
+#define misc_hausdorff_dist_normed(...) fnft__misc_hausdorff_dist_normed(__VA_ARGS__)
 #define misc_sech(...) fnft__misc_sech(__VA_ARGS__)
 #define misc_l2norm2(...) fnft__misc_l2norm2(__VA_ARGS__)
 #define misc_filter(...) fnft__misc_filter(__VA_ARGS__)
diff --git a/matlab/mex_fnft_kdvv_inverse.c b/matlab/mex_fnft_kdvv_inverse.c
new file mode 100644
index 0000000..c52187e
--- /dev/null
+++ b/matlab/mex_fnft_kdvv_inverse.c
@@ -0,0 +1,177 @@
+/*
+* This file is part of FNFT.
+*
+* FNFT is free software; you can redistribute it and/or
+* modify it under the terms of the version 2 of the GNU General
+* Public License as published by the Free Software Foundation.
+*
+* FNFT is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with this program. If not, see <http://www.gnu.org/licenses/>.
+*
+* Contributors:
+* Sander Wahls (KIT) 2026.
+*/
+
+#include <string.h>
+#include <stdio.h>
+#include "mex.h"
+#ifndef SKIP_MATRIX_H
+#include "matrix.h"
+#endif
+#include "fnft_kdvv_inverse.h"
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+    FNFT_UINT D;
+    FNFT_COMPLEX * q = NULL;
+    FNFT_REAL * T;
+    FNFT_UINT M, K;
+    FNFT_COMPLEX * contspec = NULL;
+    FNFT_COMPLEX * bound_states = NULL;
+    FNFT_COMPLEX * norming_constants = NULL;
+    FNFT_REAL * XI;
+    FNFT_UINT i;
+    double *re, *im;
+    char msg[128]; // buffer for error messages
+    FNFT_INT ret_code;
+    FNFT_UINT k;
+
+    if (nlhs < 1)
+        return;
+
+    /* Check types and dimensions of the first seven inputs: contspec, XI,
+       bound_states, norming_constants, D, T */
+
+    if ( nrhs < 6 )
+        mexErrMsgTxt("At least seven inputs expected.");
+    if ( !mxIsEmpty(prhs[0]) && (!mxIsDouble(prhs[0]) || !mxIsComplex(prhs[0]) || mxGetM(prhs[0]) != 1) )
+        mexErrMsgTxt("First input contspec should be a complex row vector (double precision) or []. Try passing complex(double(contspec(:)')).");
+    if ( mxIsComplex(prhs[1]) || !mxIsDouble(prhs[1]) || mxGetM(prhs[1]) != 1 || mxGetN(prhs[1]) != 2 )
+        mexErrMsgTxt("Second input XI should be a real 1x2 vector (double precision).");
+    if ( !mxIsEmpty(prhs[2]) && (!mxIsDouble(prhs[2]) || !mxIsComplex(prhs[2]) || mxGetM(prhs[2]) != 1) )
+        mexErrMsgTxt("Third input bound_states should be a complex row vector (double precision) or []. Try passing complex(double(bound_states(:)')).");
+    if ( !mxIsEmpty(prhs[3]) && (!mxIsDouble(prhs[3]) || !mxIsComplex(prhs[3]) || mxGetM(prhs[3]) != 1) )
+        mexErrMsgTxt("Fourth input norming_constants should be a complex row vector (double precision) or []. Try passing complex(double(norming_constants(:)')).");
+    if ( mxIsComplex(prhs[4]) || !mxIsDouble(prhs[4]) || mxGetNumberOfElements(prhs[4]) != 1 )
+        mexErrMsgTxt("Fifth input D should be a real scalar (double precision).");
+    if ( mxIsComplex(prhs[5]) || !mxIsDouble(prhs[5]) || mxGetM(prhs[5]) != 1 || mxGetN(prhs[5]) != 2 )
+        mexErrMsgTxt("Sixth input T should be a real 1x2 vector (double precision).");
+
+    M = mxGetNumberOfElements(prhs[0]);
+    K = mxGetNumberOfElements(prhs[2]);
+    T = mxGetPr(prhs[5]);
+    D = (unsigned int)mxGetScalar(prhs[4]);
+    XI = mxGetPr(prhs[1]);
+
+    /* Check values of first four inputs */
+
+    if ( K != mxGetNumberOfElements(prhs[3]) )
+        mexErrMsgTxt("bound_states and norming_constants should have the same lengths.");
+    if ( T[0] >= T[1] )
+        mexErrMsgTxt("T(1) >= T(2).");
+    if ( XI[0] >= XI[1] )
+        mexErrMsgTxt("XI(1) >= XI(2).");
+    if ( D<2 )
+        mexErrMsgTxt("D < 2.");
+
+    /* Redirect FNFT error messages and warnings to Matlabs command window */
+
+    fnft_errwarn_setprintf(mexPrintf);
+
+    /* Check remaining inputs, if any */
+
+    for (k=7; k<(FNFT_UINT)nrhs; k++) {
+
+        /* Check if current input is a string as desired and convert it */
+        if ( !mxIsChar(prhs[k]) ) {
+            snprintf(msg, sizeof msg, "%uth input should be a string.",
+                     (unsigned int)(k+1));
+            goto on_error;
+        }
+        char *str = mxArrayToString(prhs[k]);
+        if ( str == NULL ) {
+            snprintf(msg, sizeof msg, "Out of memory.");
+            goto on_error;
+        }
+
+        /* Try to interpret value of string input */
+        if ( strcmp(str, "quiet") == 0 ) {
+
+            fnft_errwarn_setprintf(NULL);
+
+        } else {
+            snprintf(msg, sizeof msg, "%uth input has invalid value.",
+                 (unsigned int)(k+1));
+            goto on_error;
+        }
+    }
+
+    /* Allocate memory */
+
+    q = mxMalloc(D * sizeof(FNFT_COMPLEX));
+    if (M>0)
+        contspec = mxMalloc(M * sizeof(FNFT_COMPLEX));
+    if (K>0) {
+        bound_states = mxMalloc(K * sizeof(FNFT_COMPLEX));
+        norming_constants = mxMalloc(K * sizeof(FNFT_COMPLEX));
+    }
+    if ( q == NULL || (M>0 && contspec == NULL) || (K>0 && bound_states == NULL)
+        || (K>0 && norming_constants == NULL) ) {
+        snprintf(msg, sizeof msg, "Out of memory.");
+        goto on_error;
+    }
+
+    /* Convert inputs */
+
+    re = mxGetPr(prhs[0]);
+    im = mxGetPi(prhs[0]);
+    for (i=0; i<M; i++)
+        contspec[i] = re[i] + I*im[i];
+    re = mxGetPr(prhs[2]);
+    im = mxGetPi(prhs[2]);
+    for (i=0; i<K; i++)
+        bound_states[i] = re[i] + I*im[i];
+    re = mxGetPr(prhs[3]);
+    im = mxGetPi(prhs[3]);
+    for (i=0; i<K; i++)
+        norming_constants[i] = re[i] + I*im[i];
+
+    /* Call the C routine */
+
+    ret_code = fnft_kdvv_inverse(M, contspec, XI, K, bound_states,
+                                 norming_constants, D, q, T, NULL);
+    if (ret_code != FNFT_SUCCESS) {
+        snprintf(msg, sizeof msg, "fnft_kdvv_inverse failed (error code %i).",
+                ret_code);
+        goto on_error;
+    }
+
+    /* Allocate memory for the output */
+
+    plhs[0] = mxCreateDoubleMatrix(1, D, mxCOMPLEX);
+
+    /* Allocate memory for outputs and convert results */
+
+    re = mxGetPr(plhs[0]);
+    im = mxGetPi(plhs[0]);
+    for (i=0; i<D; i++) {
+        re[i] = FNFT_CREAL(q[i]);
+        im[i] = FNFT_CIMAG(q[i]);
+    }
+
+    /* Free memory that is no longer needed */
+
+    mxFree(q);
+    mxFree(contspec);
+    mxFree(bound_states);
+    mxFree(norming_constants);
+    return;
+
+on_error:
+    mexErrMsgTxt(msg);
+}
diff --git a/matlab/mex_fnft_kdvv_inverse.m b/matlab/mex_fnft_kdvv_inverse.m
new file mode 100644
index 0000000..29245bc
--- /dev/null
+++ b/matlab/mex_fnft_kdvv_inverse.m
@@ -0,0 +1,51 @@
+% MEX_FNFT_KDVV_INVERSE Fast inverse nonlinear Fourier transform for the
+% Korteweg-de  Vries equation with vanishing boundaries.
+%
+%   q = MEX_FNFT_KDVVV_INVERSE(contspec, XI, bound_states, ...
+%                             norming_constants, D, T);
+%
+% DESCRIPTION
+%   Provides an interface to the C routine fnft_kdvv_inverse.
+%
+% INPUTS
+%   contspec        Complex row vector of length M>=D, contains the samples
+%                   of the reflection coefficient, the b-scattering 
+%                   coefficient or the inverse Fourier transform of the 
+%                   b-scattering coefficient on an equidistant grid.
+%                   Pass [] if the continuous spectrum is zero 
+%                   (i.e., a multi-soliton is desired)
+%                   Note: contspec functionality currently out of function 
+%                   (state 04/2026)! 
+%                   It will be implicitly assumed to \f$ 0 \f$ for all \f$ \xi \f$.
+%   XI              Real 1x2 vector, contains the position of the first and the last
+%                   sample of the continuous spectrum.
+%                   Note: contspec functionality currently out of function 
+%                   (state 04/2026)! 
+%                   It does not matter which values are chosen at the moment.
+%   bound_states    Complex row vector, contains the desired bound states.
+%   norming_constants Complex row vector, same length as bound_states.
+%                   Contains the corresponding norming constants.
+%   D               Real scalar, number of time domain samples.
+%   T               Real 1x2 vector, contains the location of the first and
+%                   the last sample in q
+%
+% OUTPUTS
+%   q               Complex row vector of length D
+
+% This file is part of FNFT.
+%
+% FNFT is free software; you can redistribute it and/or
+% modify it under the terms of the version 2 of the GNU General
+% Public License as published by the Free Software Foundation.
+%
+% FNFT is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+% Contributors:
+% Sander Wahls (KIT) 2026.
+% Fabian Fischer (Hiwi KIT) 2026.
\ No newline at end of file
diff --git a/src/fnft_kdvv_inverse.c b/src/fnft_kdvv_inverse.c
new file mode 100644
index 0000000..7c18024
--- /dev/null
+++ b/src/fnft_kdvv_inverse.c
@@ -0,0 +1,774 @@
+/*
+ * This file is part of FNFT.
+ *
+ * FNFT is free software; you can redistribute it and/or
+ * modify it under the terms of the version 2 of the GNU General
+ * Public License as published by the Free Software Foundation.
+ *
+ * FNFT is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ * Contributors:
+ * Sander Wahls (KIT) 2026
+ * Fabian Fischer (Hiwi KIT) 2026
+ * 
+ * Following Algorithms are based on Crum Transformation Matlab project 
+ * of P.J. Prins 2020
+ */
+
+
+#define FNFT_ENABLE_SHORT_NAMES
+
+#define THRESHOLD_EQUALITY_BOUND_STATES 1e-15
+
+#include "fnft_kdvv_inverse.h"
+
+static INT find_first_positive_value(
+    COMPLEX const * const array,
+    const UINT length,
+    UINT * const i_pos_ptr)
+{
+    INT ret_code = SUCCESS;
+    UINT is_positive_value_found = 0;
+    
+    for (UINT i = 0; i < length; i++){
+        if (CREAL(array[i]) > 0){
+            *i_pos_ptr = i;
+            is_positive_value_found = 1;
+            break;
+        }
+    }
+
+    if (is_positive_value_found == 0){
+        ret_code = FNFT_EC_INVALID_ARGUMENT;
+    }
+    else {
+        ret_code = SUCCESS;
+    }
+
+    return ret_code;
+}
+
+static INT one_solition_crum_transformation(
+    COMPLEX const k1,
+    COMPLEX const * const th1,
+    COMPLEX const * const th2,
+    const UINT D,
+    const UINT i_pos,
+    COMPLEX const * const t_grid,
+    COMPLEX * const w_inv,
+    COMPLEX * const prefactor,
+    COMPLEX * const M_min1_11)
+{
+    INT ret_code = SUCCESS;
+
+    // Calculation for positive x
+    for (UINT i = i_pos; i<D; i++){
+        w_inv[i] = 1/(th1[i] + th2[i] * CEXP(2*k1*t_grid[i]));
+        prefactor[i] = 2 * k1*k1 * th1[i] * th2[i] * w_inv[i]*w_inv[i];
+        M_min1_11[i] = prefactor[i] * CEXP(2*k1*t_grid[i]);
+    }    
+
+    // Calculation for negative x
+    for (UINT i = 0; i<i_pos; i++){
+        w_inv[i] = 1/(th1[i] * CEXP(-2*k1*t_grid[i]) + th2[i]);
+        prefactor[i] = 2 * k1*k1 * th1[i] * th2[i] * w_inv[i]*w_inv[i];
+        M_min1_11[i] = prefactor[i] *  CEXP(-2*k1*t_grid[i]);
+    }
+
+    return ret_code;
+}
+
+static INT two_solitions_crum_transformation(
+    COMPLEX const k1,
+    COMPLEX const k2,
+    const UINT N_pm0_jZ,
+    COMPLEX const * const pm0_jZ,
+    COMPLEX const * const theta_E1,
+    COMPLEX const * const theta_E2,
+    COMPLEX * const tmp_t1,
+    COMPLEX * const tmp_t2,
+    const UINT D,
+    const UINT i_pos,
+    COMPLEX const * const t_grid,
+    COMPLEX * const w_inv,
+    COMPLEX * const prefactor,
+    COMPLEX * const dq,
+    COMPLEX * const s)
+{
+    INT ret_code = SUCCESS;
+
+    // Defining theta vectors out of theta_E vectors, belonging to the bound states to add
+    // These vectors are used in the Crum-Transformation step
+    // thetas according to 1. eigenvalue
+    COMPLEX const * const th11 = &theta_E1[0];
+    COMPLEX const * const th12 = &theta_E2[0];
+
+    // thetas according to 2. eigenvalue
+    COMPLEX const * const th21 = &theta_E1[D];
+    COMPLEX const * const th22 = &theta_E2[D];
+
+    // for positive x
+    for (UINT i=i_pos; i<D; i++){
+        tmp_t1[i] = th11[i] + th12[i] * CEXP(-2*k1*t_grid[i]);
+        tmp_t2[i] = th21[i] + th22[i] * CEXP(-2*k2*t_grid[i]);
+        w_inv[i] = 1/(tmp_t1[i] * k2 * (th21[i] - th22[i] * CEXP(-2*k2*t_grid[i])) - 
+                    tmp_t2[i] * k1 * (th11[i] - th12[i] * CEXP(-2*k1*t_grid[i])));
+        prefactor[i] = 2 * (k2*k2 - k1*k1) * w_inv[i]*w_inv[i];
+
+        for (UINT j=0; j<N_pm0_jZ; j++) {
+            s[j*D+i] =  tmp_t2[i]*tmp_t2[i] * k1*k1 * th11[i] * th12[i] * CEXP(-2*(k1-pm0_jZ[j])*t_grid[i]) -
+                        tmp_t1[i]*tmp_t1[i] * k2*k2 * th21[i] * th22[i] * CEXP(-2*(k2-pm0_jZ[j])*t_grid[i]);
+        }
+
+        dq[i] = prefactor[i] * s[(N_pm0_jZ-1)*D+i];
+    }
+
+    // for negative x
+    for (UINT i=0; i<i_pos; i++){
+        tmp_t1[i] = th11[i] * CEXP(2*k1*t_grid[i]) + th12[i];
+        tmp_t2[i] = th21[i] * CEXP(2*k2*t_grid[i]) + th22[i];
+        w_inv[i] = 1/(tmp_t1[i] * k2 * (th21[i] * CEXP(2*k2*t_grid[i]) - th22[i]) - 
+                    tmp_t2[i] * k1 * (th11[i] * CEXP(2*k1*t_grid[i]) - th12[i]));
+        prefactor[i] = 2 * (k2*k2 - k1*k1) * w_inv[i]*w_inv[i];
+
+        for (UINT j=0; j<N_pm0_jZ; j++) {
+            s[j*D+i] =  tmp_t2[i]*tmp_t2[i] * k1*k1 * th11[i] * th12[i] * CEXP(2*(k1+pm0_jZ[j])*t_grid[i]) -
+                        tmp_t1[i]*tmp_t1[i] * k2*k2 * th21[i] * th22[i] * CEXP(2*(k2+pm0_jZ[j])*t_grid[i]);
+        }
+
+        dq[i] = prefactor[i] * s[(N_pm0_jZ-1)*D+i];
+    }
+
+    return ret_code;
+}
+
+
+
+static INT one_solition_update_Jost(
+    const COMPLEX k1,
+    const UINT N_jZ,
+    COMPLEX const * const bound_states_left,
+    COMPLEX * const theta_E1,
+    COMPLEX * const theta_E2,
+    const UINT D,
+    const UINT i_pos,
+    COMPLEX const * const t_grid,
+    COMPLEX const * const w_inv,
+    COMPLEX const * const prefactor,
+    COMPLEX const * const M_min1_11)
+{
+    INT ret_code = SUCCESS;
+    
+    COMPLEX * C_E = NULL;
+    COMPLEX * prefactor_C_E = NULL;
+
+    // Check, if Jost has to be updated. Only when there are bound states left, that are not added yet
+    UINT const is_Jost_to_update = N_jZ != 0;
+
+    // Defining theta vectors out of theta_E vectors, belonging to the bound state to add
+    COMPLEX const * const th1 = &theta_E1[0];
+    COMPLEX const * const th2 = &theta_E2[0];
+
+    // theta vectors according to the bound states left, which corresponds to the Jost solutions
+    COMPLEX * const th1_jZ = &theta_E1[D];
+    COMPLEX * const th2_jZ = &theta_E2[D];
+
+    // Create C_E Matrix (for Updating Jost Solution)
+    UINT n_C_E = 2*2*D*N_jZ; 
+    C_E = malloc(n_C_E*sizeof(COMPLEX));
+    CHECK_NOMEM(C_E, ret_code, leave_fun);
+
+    for (UINT i=0; i<n_C_E; i++){
+        C_E[i] = 0;
+    }
+
+    COMPLEX * C_E_1_1 = &C_E[0*D*N_jZ];
+    COMPLEX * C_E_1_2 = &C_E[1*D*N_jZ];
+    COMPLEX * C_E_2_1 = &C_E[2*D*N_jZ];
+    COMPLEX * C_E_2_2 = &C_E[3*D*N_jZ];
+
+    for (UINT i=0; i<D; i++){
+        for (UINT j=0; j<N_jZ; j++){
+            C_E_1_1[j*D+i] = -bound_states_left[j];
+            C_E_2_2[j*D+i] = bound_states_left[j];
+        }
+    }
+
+    // Compute a new prefactor for C_E
+    prefactor_C_E = malloc(D * sizeof(COMPLEX));
+    CHECK_NOMEM(prefactor_C_E, ret_code, leave_fun);
+
+    // Calculation for positive x
+    for (UINT i=i_pos; i<D && is_Jost_to_update; i++){
+        prefactor_C_E[i] = k1 * (th1[i] - th2[i] * CEXP(2*k1*t_grid[i])) * w_inv[i];
+
+        for (UINT j=0; j<N_jZ; j++){
+            C_E_1_1[j*D+i] = C_E_1_1[j*D+i] + prefactor_C_E[i] + M_min1_11[i] * 1/bound_states_left[j];
+            C_E_1_2[j*D+i] = C_E_1_2[j*D+i] + prefactor[i] * CEXP(2*(k1+bound_states_left[j])*t_grid[i])*1/bound_states_left[j];
+            C_E_2_1[j*D+i] = C_E_2_1[j*D+i] - prefactor[i] * CEXP(2*(k1-bound_states_left[j])*t_grid[i])*1/bound_states_left[j];
+            C_E_2_2[j*D+i] = C_E_2_2[j*D+i] + prefactor_C_E[i] - M_min1_11[i] * 1/bound_states_left[j];
+        }
+    }
+    
+    // Calculation for negative x
+    for (UINT i=0; i<i_pos && is_Jost_to_update; i++){
+        prefactor_C_E[i] = k1 * (th1[i] * CEXP(-2*k1*t_grid[i]) - th2[i]) * w_inv[i];
+        
+        for (UINT j=0; j<N_jZ; j++){
+            C_E_1_1[j*D+i] = C_E_1_1[j*D+i] + prefactor_C_E[i] + M_min1_11[i] * 1/bound_states_left[j];
+            C_E_1_2[j*D+i] = C_E_1_2[j*D+i] + prefactor[i] * CEXP(-2*(k1-bound_states_left[j])*t_grid[i])*1/bound_states_left[j];
+            C_E_2_1[j*D+i] = C_E_2_1[j*D+i] - prefactor[i] * CEXP(-2*(k1+bound_states_left[j])*t_grid[i])*1/bound_states_left[j];
+            C_E_2_2[j*D+i] = C_E_2_2[j*D+i] + prefactor_C_E[i] - M_min1_11[i] * 1/bound_states_left[j];
+        }
+    }
+    
+    // Map Jost solution
+    COMPLEX tmp_th1;
+    COMPLEX tmp_th2;
+
+    for (UINT i=0; i<D && is_Jost_to_update; i++){
+        for (UINT j=0; j<N_jZ; j++){
+            tmp_th1 = th1_jZ[j*D+i];
+            tmp_th2 = th2_jZ[j*D+i];
+
+            th1_jZ[j*D+i] = C_E_1_1[j*D+i] * tmp_th1 + C_E_1_2[j*D+i] * tmp_th2;
+            th2_jZ[j*D+i] = C_E_2_1[j*D+i] * tmp_th1 + C_E_2_2[j*D+i] * tmp_th2;
+        }
+    }
+
+leave_fun:
+    free(C_E);
+    free(prefactor_C_E);
+
+    return ret_code;
+}
+
+static INT two_solitions_update_Jost(
+    const COMPLEX k1,
+    const COMPLEX k2,
+    const UINT N_jZ,
+    COMPLEX const * const bound_states_left,
+    COMPLEX const * const pm0_jZ,
+    COMPLEX * const theta_E1,
+    COMPLEX * const theta_E2,
+    COMPLEX const * const temp_t1,
+    COMPLEX const * const temp_t2,
+    const UINT D,
+    const UINT i_pos,
+    COMPLEX const * const t_grid,
+    COMPLEX const * const w_inv,
+    COMPLEX const * const prefactor,
+    COMPLEX const * const s)
+{
+    INT ret_code = SUCCESS;
+    COMPLEX * C_E = NULL;
+
+    // Check, if Jost has to be updated. Only when there are bound states left, that are not added yet
+    UINT const is_Jost_to_update = N_jZ != 0;
+
+    // Defining theta vectors out of theta_E vectors, belonging to the bound states to add
+    // These vectors are used in the Crum-Transformation step
+    // thetas according to 1. eigenvalue
+    COMPLEX const * const th11 = &theta_E1[0];
+    COMPLEX const * const th12 = &theta_E2[0];
+
+    // thetas according to 2. eigenvalue
+    COMPLEX const * const th21 = &theta_E1[D];
+    COMPLEX const * const th22 = &theta_E2[D];
+
+    // theta vectors according to the bound states left, which corresponds to the Jost solutions
+    COMPLEX * const th1_jZ = &theta_E1[2*D];
+    COMPLEX * const th2_jZ = &theta_E2[2*D];
+
+    // Create C_E Matrix (for Updating Jost Solution)
+    UINT n_C_E = 2*2*D*N_jZ; 
+    C_E = malloc(n_C_E*sizeof(COMPLEX));
+    CHECK_NOMEM(C_E, ret_code, leave_fun);
+
+    for (UINT i=0; i<n_C_E; i++){
+        C_E[i] = 0;
+    }
+
+    COMPLEX * C_E_1_1 = &C_E[0*D*N_jZ];
+    COMPLEX * C_E_1_2 = &C_E[1*D*N_jZ];
+    COMPLEX * C_E_2_1 = &C_E[2*D*N_jZ];
+    COMPLEX * C_E_2_2 = &C_E[3*D*N_jZ];
+
+    for (UINT i=0; i<D; i++){
+        for (UINT j=0; j<N_jZ; j++){
+            C_E_1_1[j*D+i] = bound_states_left[j]*bound_states_left[j];
+            C_E_2_2[j*D+i] = bound_states_left[j]*bound_states_left[j];
+        }
+    }
+
+    COMPLEX prefactor1;
+    COMPLEX M_0_11;
+    COMPLEX prefactor2;
+    COMPLEX mpjZinv_0;
+    COMPLEX mpjZinv_1;
+    COMPLEX mpjZinv_2;
+    COMPLEX M_min1_11_jZinv;
+
+    // Calculation for positive x
+    for (UINT i=i_pos; i<D && is_Jost_to_update; i++){
+        prefactor1 = (k2*k2 - k1*k1) * temp_t1[i] * temp_t2[i] * w_inv[i];
+        M_0_11 = -(k1*k1 + k2*k2)/2 + prefactor[i] * (k2*k2 - k1*k1)/4 * temp_t1[i]*temp_t1[i] * temp_t2[i]*temp_t2[i];
+        prefactor2 = prefactor[i] * k1 * k2;
+        
+        mpjZinv_0 = (k2 * th21[i] * th22[i] * (th11[i]*th11[i] - th12[i]*th12[i] * 
+                    CEXP(-4*k1*t_grid[i])) * CEXP(-2*(k2) * t_grid[i]) + 
+                    -k1 * th11[i] * th12[i] * (th21[i]*th21[i] - th22[i]*th22[i] * 
+                    CEXP(-4*k2*t_grid[i])) * CEXP(-2*(k1) * t_grid[i]));
+
+        for (UINT j=0; j<N_jZ; j++){
+            mpjZinv_1 = (k2 * th21[i] * th22[i] * (th11[i]*th11[i] - th12[i]*th12[i] * 
+                        CEXP(-4*k1*t_grid[i])) * CEXP(-2*(k2-pm0_jZ[j]) * t_grid[i]) + 
+                        -k1 * th11[i] * th12[i] * (th21[i]*th21[i] - th22[i]*th22[i] * 
+                        CEXP(-4*k2*t_grid[i])) * CEXP(-2*(k1-pm0_jZ[j]) * t_grid[i])  
+                    ) * 1/bound_states_left[j];
+
+            mpjZinv_2 = (k2 * th21[i] * th22[i] * (th11[i]*th11[i] - th12[i]*th12[i] * 
+                        CEXP(-4*k1*t_grid[i])) * CEXP(-2*(k2-pm0_jZ[j+N_jZ]) * t_grid[i]) + 
+                        -k1 * th11[i] * th12[i] * (th21[i]*th21[i] - th22[i]*th22[i] * 
+                        CEXP(-4*k2*t_grid[i])) * CEXP(-2*(k1-pm0_jZ[j+N_jZ]) * t_grid[i])  
+                    ) * 1/bound_states_left[j];
+
+            M_min1_11_jZinv = prefactor2 * mpjZinv_0 * 1/bound_states_left[j];
+
+            C_E_1_1[j*D+i] = C_E_1_1[j*D+i] + prefactor1 * bound_states_left[j] + M_0_11 + M_min1_11_jZinv;
+            C_E_1_2[j*D+i] = C_E_1_2[j*D+i] + prefactor[i] * s[j*D+i] + prefactor2 * mpjZinv_1;
+            C_E_2_1[j*D+i] = C_E_2_1[j*D+i] + prefactor[i] * s[(j+N_jZ)*D+i] - prefactor2 * mpjZinv_2;
+            C_E_2_2[j*D+i] = C_E_2_2[j*D+i] - prefactor1 * bound_states_left[j] + M_0_11 - M_min1_11_jZinv;
+        }
+    }
+
+    // Calculation for negative x
+    for (UINT i=0; i<i_pos && is_Jost_to_update; i++){
+        prefactor1 = (k2*k2 - k1*k1) * temp_t1[i] * temp_t2[i] * w_inv[i];
+        M_0_11 = -(k1*k1 + k2*k2)/2 + prefactor[i] * (k2*k2 - k1*k1)/4 * temp_t1[i]*temp_t1[i] * temp_t2[i]*temp_t2[i];
+        prefactor2 = prefactor[i] * k1 * k2;
+
+        mpjZinv_0 = (k2 * th21[i] * th22[i] * (th11[i]*th11[i] * CEXP(4*k1*t_grid[i]) - th12[i]*th12[i]) *               
+                    CEXP(2*(k2) * t_grid[i]) +
+                    -k1 * th11[i] * th12[i] * (th21[i]*th21[i] * CEXP(4*k2*t_grid[i]) - th22[i]*th22[i]) *
+                    CEXP(2*(k1) * t_grid[i]) );
+
+        for (UINT j=0; j<N_jZ; j++){
+            mpjZinv_1 = (k2 * th21[i] * th22[i] * (th11[i]*th11[i] * CEXP(4*k1*t_grid[i]) - 
+                        th12[i]*th12[i]) * CEXP(2*(k2+pm0_jZ[j]) * t_grid[i]) +
+                        -k1 * th11[i] * th12[i] * (th21[i]*th21[i] * CEXP(4*k2*t_grid[i]) - 
+                        th22[i]*th22[i]) * CEXP(2*(k1+pm0_jZ[j]) * t_grid[i])
+                    ) * 1/bound_states_left[j];
+
+            mpjZinv_2 = (k2 * th21[i] * th22[i] * (th11[i]*th11[i] * CEXP(4*k1*t_grid[i]) - 
+                        th12[i]*th12[i]) * CEXP(2*(k2+pm0_jZ[j+N_jZ]) * t_grid[i]) +
+                        -k1 * th11[i] * th12[i] * (th21[i]*th21[i] * CEXP(4*k2*t_grid[i]) - 
+                        th22[i]*th22[i]) * CEXP(2*(k1+pm0_jZ[j+N_jZ]) * t_grid[i])
+                    ) * 1/bound_states_left[j];
+
+            M_min1_11_jZinv = prefactor2 * mpjZinv_0 * 1/bound_states_left[j];
+
+            C_E_1_1[j*D+i] = C_E_1_1[j*D+i] + prefactor1 * bound_states_left[j] + M_0_11 + M_min1_11_jZinv;
+            C_E_1_2[j*D+i] = C_E_1_2[j*D+i] + prefactor[i] * s[j*D+i] + prefactor2 * mpjZinv_1;
+            C_E_2_1[j*D+i] = C_E_2_1[j*D+i] + prefactor[i] * s[(j+N_jZ)*D+i] - prefactor2 * mpjZinv_2;
+            C_E_2_2[j*D+i] = C_E_2_2[j*D+i] - prefactor1 * bound_states_left[j] + M_0_11 - M_min1_11_jZinv;
+        }
+    }
+
+    // Map Jost solution
+    COMPLEX tmp_th1;
+    COMPLEX tmp_th2;
+
+    for (UINT i=0; i<D && is_Jost_to_update; i++){
+        for (UINT j=0; j<N_jZ; j++){
+            tmp_th1 = th1_jZ[j*D+i];
+            tmp_th2 = th2_jZ[j*D+i];
+
+            th1_jZ[j*D+i] = C_E_1_1[j*D+i] * tmp_th1 + C_E_1_2[j*D+i] * tmp_th2;
+            th2_jZ[j*D+i] = C_E_2_1[j*D+i] * tmp_th1 + C_E_2_2[j*D+i] * tmp_th2;
+        }
+    }
+
+leave_fun:
+    free(C_E);
+
+    return ret_code;
+}
+
+
+static INT add_one_soliton(
+    COMPLEX const * const bound_states,
+    UINT const N_bound_states,
+    COMPLEX * const theta_E1,
+    COMPLEX * const theta_E2,
+    COMPLEX const * const t_grid,
+    const UINT D,
+    COMPLEX * const q)
+{
+    // Initialize variables and arrays
+    INT ret_code = SUCCESS;
+
+    COMPLEX * M_min1_11 = NULL;
+    COMPLEX * w_inv = NULL;
+    COMPLEX * prefactor = NULL;
+
+    M_min1_11 = malloc(D * sizeof(COMPLEX));
+    CHECK_NOMEM(M_min1_11, ret_code, leave_fun);
+    w_inv = malloc(D * sizeof(COMPLEX));
+    CHECK_NOMEM(w_inv, ret_code, leave_fun);
+    prefactor = malloc(D * sizeof(COMPLEX));
+    CHECK_NOMEM(prefactor, ret_code, leave_fun);
+
+    // Determine where x>0
+    UINT i_pos;
+    ret_code = find_first_positive_value(t_grid, D, &i_pos);
+    CHECK_RETCODE(ret_code, leave_fun);
+
+    // -- Crum Transform --
+    COMPLEX bound_state_added_here = bound_states[0];
+    ret_code = one_solition_crum_transformation(    bound_state_added_here,
+                                                    &theta_E1[0],           // thetas belonging to the bound state to add
+                                                    &theta_E2[0],           // thetas belonging to the bound state to add
+                                                    D,
+                                                    i_pos,
+                                                    t_grid,
+                                                    w_inv,
+                                                    prefactor,
+                                                    M_min1_11);
+    CHECK_RETCODE(ret_code, leave_fun);                                                
+    
+    // Update output
+    for (UINT i=0; i<D; i++){
+        q[i] = -q[i] + 4 * M_min1_11[i];      
+    }
+    
+    // Update Jost Solution only when there are further bound states left, which are not added yet
+    UINT N_bound_states_left = N_bound_states - 1; 
+    if (N_bound_states_left) {
+        COMPLEX const * const bound_states_left = &bound_states[1];
+        ret_code = one_solition_update_Jost(    bound_state_added_here,
+                                                N_bound_states_left,        
+                                                bound_states_left,
+                                                theta_E1,
+                                                theta_E2,
+                                                D,
+                                                i_pos,
+                                                t_grid,
+                                                w_inv,
+                                                prefactor,
+                                                M_min1_11);
+        CHECK_RETCODE(ret_code, leave_fun);                                            
+    }
+
+leave_fun:
+    free(M_min1_11);
+    free(w_inv);
+    free(prefactor);
+
+    return ret_code;
+}
+
+
+static INT add_two_solitons(
+    COMPLEX const * const bound_states,
+    UINT const N_bound_states,
+    COMPLEX * const theta_E1,
+    COMPLEX * const theta_E2,
+    COMPLEX const * const t_grid,
+    const UINT D,
+    COMPLEX * const q)
+{
+    INT ret_code = SUCCESS;
+
+    COMPLEX * w_inv = NULL;
+    COMPLEX * prefactor = NULL;
+    COMPLEX * dq = NULL;
+    COMPLEX * pm0_jZ = NULL;
+    COMPLEX * s = NULL;
+    COMPLEX * temp_t1 = NULL;
+    COMPLEX * temp_t2 = NULL;
+
+    w_inv = malloc(D * sizeof(COMPLEX));
+    CHECK_NOMEM(w_inv, ret_code, leave_fun);
+    prefactor = malloc(D * sizeof(COMPLEX));
+    CHECK_NOMEM(prefactor, ret_code, leave_fun);
+    dq = malloc(D * sizeof(COMPLEX));
+    CHECK_NOMEM(dq, ret_code, leave_fun);
+
+    // For the bound states, which are not 
+    const UINT N_bound_states_left = N_bound_states - 2;
+    UINT const is_Jost_to_update = N_bound_states_left != 0;
+
+    COMPLEX const * bound_states_left = NULL;
+
+    if (is_Jost_to_update) {
+        bound_states_left = &bound_states[2];
+    }
+
+    // Initialize pm0_jZ
+    UINT N_pm0_jZ = 2 * N_bound_states_left + 1;
+    pm0_jZ = malloc((2*N_bound_states_left+1) * sizeof(COMPLEX));
+    CHECK_NOMEM(pm0_jZ, ret_code, leave_fun);
+    
+    for (UINT i=0; i<N_bound_states_left && is_Jost_to_update; i++){
+        pm0_jZ[i] = bound_states_left[i];
+        pm0_jZ[i+N_bound_states_left] = -bound_states_left[i];
+    }
+    
+    pm0_jZ[2*N_bound_states_left] = 0;
+
+    // dimension of s is dependend on number of bound states left
+    s = malloc(D * N_pm0_jZ * sizeof(COMPLEX));
+    CHECK_NOMEM(s, ret_code, leave_fun);
+
+    // Intermediate variables
+    temp_t1 = malloc(D * sizeof(COMPLEX));
+    CHECK_NOMEM(temp_t1, ret_code, leave_fun);
+    temp_t2 = malloc(D * sizeof(COMPLEX));
+    CHECK_NOMEM(temp_t2, ret_code, leave_fun);
+    
+    // Change Sign of the bound states to add
+    COMPLEX k1 = -bound_states[0];
+    COMPLEX k2 = -bound_states[1];
+
+    // Determine where x>0
+    UINT i_pos;
+    ret_code = find_first_positive_value(t_grid, D, &i_pos);
+    CHECK_RETCODE(ret_code, leave_fun);
+
+    // -- Crum-Transformation --
+    ret_code = two_solitions_crum_transformation(   k1,
+                                                    k2,
+                                                    N_pm0_jZ,
+                                                    pm0_jZ,
+                                                    theta_E1,
+                                                    theta_E2,
+                                                    temp_t1,
+                                                    temp_t2,
+                                                    D,
+                                                    i_pos,
+                                                    t_grid,
+                                                    w_inv,
+                                                    prefactor,
+                                                    dq,
+                                                    s   );
+    CHECK_RETCODE(ret_code, leave_fun);
+
+
+    // Update output
+    for (UINT i=0; i<D; i++){
+        q[i] = q[i] + 4 * dq[i];      
+    }
+
+    // -- Update Jost Solution --
+    if (is_Jost_to_update) {
+        ret_code = two_solitions_update_Jost(   k1,
+                                                k2,
+                                                N_bound_states_left,
+                                                bound_states_left,
+                                                pm0_jZ,
+                                                theta_E1,
+                                                theta_E2,
+                                                temp_t1,
+                                                temp_t2,
+                                                D,
+                                                i_pos,
+                                                t_grid,
+                                                w_inv,
+                                                prefactor,
+                                                s   );
+        CHECK_RETCODE(ret_code, leave_fun);
+    }
+
+leave_fun:
+    free(w_inv);
+    free(prefactor);
+    free(dq);
+    free(pm0_jZ);
+    free(s);
+    free(temp_t1);
+    free(temp_t2);
+    
+    return ret_code;
+}
+
+
+
+INT fnft_kdvv_inverse(
+    UINT const M,
+    COMPLEX * const contspec,
+    REAL const * const XI,
+    UINT const K,
+    COMPLEX const * const bound_states,
+    COMPLEX const * const norming_constants,
+    UINT const D,
+    COMPLEX * const q,
+    REAL const * const T,
+    void * opts_ptr)
+{  
+   // checks for valid input ---------------------------------------------------------------------------------
+   
+   if (T[1] <= T[0]) {
+       return E_INVALID_ARGUMENT_MSG(T,Position of the first sample is bigger than of the last sample (i.e. T[1] <= T[0])!);
+    }
+    
+    for (UINT i=0; i<K; i++) {
+        if (CREAL(bound_states[i]) != 0){
+            return E_INVALID_ARGUMENT_MSG(bound_states,At least one bound state has a nonzero real part!);
+        }
+        
+        if (CIMAG(bound_states[i]) < 0){
+            return E_INVALID_ARGUMENT_MSG(bound_states,At least one bound state is negative imaginary!);
+        }
+        
+        if (norming_constants[i] == 0){
+            return E_INVALID_ARGUMENT_MSG(norming_constants,At least one norming constant is zero!);
+        }
+    }
+    
+    for (UINT i=0; i<K-1; i++) {
+        if (CABS(bound_states[i]) < CABS(bound_states[i+1])){
+            return E_INVALID_ARGUMENT_MSG(bound_states,The bound states have not a descend order!);
+        }
+        
+        if (CREAL(norming_constants[i]*norming_constants[i+1]) > 0){
+            return E_INVALID_ARGUMENT_MSG(norming_constants,The signs of the norming constants does not alternate!);
+        }
+
+        for (UINT j=i+1; j<K-1; j++) {
+            if (CABS(bound_states[i] - bound_states[j]) < THRESHOLD_EQUALITY_BOUND_STATES){
+                return E_INVALID_ARGUMENT_MSG(bound_states,At least two bound states are almost equal!);
+            }
+        }
+    }
+    
+    // first checking if correct order of bound states exist, than checking if the first norming constant is positive!
+    if (CREAL(norming_constants[0]) < 0){
+        return E_INVALID_ARGUMENT_MSG(norming_constants,The norming constant of the biggest bound state is negative!);
+    }
+
+    // ----------------------------------------------------------------------------------------------------------
+    
+    INT ret_code = SUCCESS;
+    
+    COMPLEX * t_grid = NULL;
+    COMPLEX * bound_states_i = NULL;
+    COMPLEX * normconsts_i = NULL;
+    COMPLEX * theta_E1 = NULL;
+    COMPLEX * theta_E2 = NULL;
+
+    // Initialize q
+    for (UINT n=0; n<D; n++){
+        q[n] = 0.0;
+    }
+
+    // Initialize spatial grid with D points between T[0] and T[1]
+    t_grid = malloc(D * sizeof(COMPLEX));
+    CHECK_NOMEM(t_grid, ret_code, leave_fun);
+
+    const COMPLEX eps_t = (T[1] - T[0])/(D - 1);
+
+    for (UINT n=0; n<D; n++) {
+        t_grid[n]= T[0] + n*eps_t;
+    }
+
+    // initialize new storage for bound states and normconsts for later sorting
+    bound_states_i = malloc(K * sizeof(COMPLEX));
+    CHECK_NOMEM(bound_states_i, ret_code, leave_fun);
+
+    normconsts_i = malloc(K * sizeof(COMPLEX));
+    CHECK_NOMEM(normconsts_i, ret_code, leave_fun);
+
+    for (UINT i=0; i<K; i++){
+        bound_states_i[i] = bound_states[i]*I; // bound states should be real for further computing
+        normconsts_i[i] = norming_constants[i];
+    }
+
+    // shift sign of norm_consts if K even
+    if (K%2 == 0){
+        for (UINT i=0; i<K; i++){
+            normconsts_i[i] = -1*normconsts_i[i];
+        }
+    } 
+
+    // Declare theta_E
+    theta_E1 = malloc(D * K * sizeof(COMPLEX));
+    CHECK_NOMEM(theta_E1, ret_code, leave_fun);
+
+    theta_E2 = malloc(D * K * sizeof(COMPLEX));
+    CHECK_NOMEM(theta_E2, ret_code, leave_fun);
+
+    for (UINT i=0; i<D; i++){
+        for (UINT j=0; j<K; j++){
+            theta_E1[j*D+i] = 1;
+            theta_E2[j*D+i] = normconsts_i[j];
+        }
+    }
+
+    UINT N_rest = K;
+    UINT step_idx = 0;
+    UINT N_step;
+
+    // Add solitions
+    while (N_rest > 0){
+        // Select indices for this step and for the remaining eigenvalues to add
+        // If the number of eigenvalues left is odd, then only one solition should be added
+        if (N_rest%2==1){ N_step = 1; }
+        else { N_step = 2; }
+        
+        // Scale trajectories (magnitudes of th1 and th2 symmetric around 0)
+        COMPLEX tmp_scaling_factor;
+
+        for (UINT i=0; i<D; i++){
+            for (UINT j=step_idx; j<step_idx+N_step; j++){
+                tmp_scaling_factor = CPOW(2, -ROUND( (LOG2(CABS(theta_E1[j*D+i])) + LOG2(CABS(theta_E2[j*D+i])))/2 ));
+                theta_E1[j*D+i] = theta_E1[j*D+i] * tmp_scaling_factor;
+                theta_E2[j*D+i] = theta_E2[j*D+i] * tmp_scaling_factor;
+            }
+        }
+
+        // Crum-Transformation step
+        if (N_step==1){
+            ret_code = add_one_soliton( &bound_states_i[step_idx], 
+                                        N_rest,
+                                        &theta_E1[step_idx*D],
+                                        &theta_E2[step_idx*D],
+                                        t_grid, 
+                                        D, 
+                                        q);
+            CHECK_RETCODE(ret_code, leave_fun);
+        } else {
+            ret_code = add_two_solitons(&bound_states_i[step_idx], 
+                                        N_rest,
+                                        &theta_E1[step_idx*D],
+                                        &theta_E2[step_idx*D],
+                                        t_grid, 
+                                        D, 
+                                        q);
+            CHECK_RETCODE(ret_code, leave_fun);
+        }
+
+        N_rest = N_rest-N_step;
+        step_idx = step_idx + N_step;
+
+        // Filter out errors that result in negative values
+        for (UINT i=0; i<D; i++){
+            if (CREAL(q[i]) < 0.0){ q[i] = 0.0; }
+        }
+    }
+
+leave_fun:
+    free(t_grid);
+    free(bound_states_i);
+    free(normconsts_i);
+    free(theta_E1);
+    free(theta_E2);
+    return ret_code;
+}
+
diff --git a/src/private/fnft__kdvv_inverse_testcases.c b/src/private/fnft__kdvv_inverse_testcases.c
new file mode 100644
index 0000000..24c6fa8
--- /dev/null
+++ b/src/private/fnft__kdvv_inverse_testcases.c
@@ -0,0 +1,192 @@
+/*
+* This file is part of FNFT.  
+*                                                                  
+* FNFT is free software; you can redistribute it and/or
+* modify it under the terms of the version 2 of the GNU General
+* Public License as published by the Free Software Foundation.
+*
+* FNFT is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*                                                                      
+* You should have received a copy of the GNU General Public License
+* along with this program. If not, see <http://www.gnu.org/licenses/>.
+*
+* Contributors:
+* Fabian Fischer (Hiwi KIT) 2026
+*/
+
+#define FNFT_ENABLE_SHORT_NAMES
+
+#include <stdio.h>
+
+#include "fnft__kdvv_inverse_testcases.h"
+#include "fnft__errwarn.h"
+#include "fnft__misc.h"
+
+void kdvv_print_spectrum(   COMPLEX const * const bound_states,
+                            COMPLEX const * const normconsts,
+                            COMPLEX const * const contspec,
+                            REAL const * const XI,
+                            const UINT M,
+                            const UINT D,
+                            const UINT K)
+{
+    printf("Number of samples:\n  D = %u\n", (unsigned int)D);
+
+    FNFT_REAL eps_xi = (XI[1] - XI[0]) / (M - 1);
+    printf("Continuous spectrum:\n");
+    for (FNFT_UINT i=0; i<M; i++) {
+        FNFT_REAL xi = XI[0] + i*eps_xi;
+        printf("  continuous_spectrum(xi=%f) \t= %g + %gI\n",
+            (double)xi,
+            (double)FNFT_CREAL(contspec[i]),
+            (double)FNFT_CIMAG(contspec[i])
+        );
+    }
+
+    printf("Discrete spectrum:\n");
+    for (FNFT_UINT i=0; i<K; i++) {
+        printf("  bound state at %g + %gI with norming constant %g + %gI\n",
+            (double)FNFT_CREAL(bound_states[i]),
+            (double)FNFT_CIMAG(bound_states[i]),
+            (double)FNFT_CREAL(normconsts[i]),
+            (double)FNFT_CIMAG(normconsts[i])
+        );
+    }
+}
+
+
+
+INT kdvv_testcases_get_spectrum_of_inverse(const fnft_kdvv_params params_i,
+                    const REAL err_bnd_bound_states,
+                    const REAL err_bnd_spurious_bound_states,
+                    const REAL err_bnd_normconst,
+                    const REAL err_bnd_contspec)
+{
+    INT ret_code = SUCCESS;
+    
+    COMPLEX * q = NULL;
+    COMPLEX * contspec_r = NULL;
+    COMPLEX * bound_states_r = NULL;
+    COMPLEX * normconsts_r = NULL;
+    
+
+    q = malloc(params_i.D * sizeof(COMPLEX));
+    CHECK_NOMEM(q, ret_code, leave_fun);
+
+    ret_code = fnft_kdvv_inverse(   params_i.M, 
+                                    params_i.contspec, 
+                                    params_i.XI, 
+                                    params_i.K, 
+                                    params_i.bound_states, 
+                                    params_i.normconsts, 
+                                    params_i.D, 
+                                    q, 
+                                    params_i.T, NULL);
+    CHECK_RETCODE(ret_code, leave_fun);
+    
+    // Prepare forward fnft_kdvv
+    contspec_r = malloc(params_i.M * sizeof(COMPLEX));
+    CHECK_NOMEM(contspec_r, ret_code, leave_fun);
+
+    UINT K_r = params_i.D;
+
+    bound_states_r = malloc(K_r * sizeof(COMPLEX));
+    CHECK_NOMEM(bound_states_r, ret_code, leave_fun);
+
+    normconsts_r = malloc(K_r * sizeof(COMPLEX));
+    CHECK_NOMEM(normconsts_r, ret_code, leave_fun);
+
+    fnft_kdvv_opts_t opts = fnft_kdvv_default_opts();
+
+    ret_code = fnft_kdvv(params_i.D, q, params_i.T, params_i.M, contspec_r, params_i.XI, &K_r, bound_states_r, normconsts_r, &opts);
+    CHECK_RETCODE(ret_code, leave_fun);
+
+    #ifdef DEBUG
+        kdvv_print_spectrum(bound_states_r, normconsts_r, contspec_r, params_i.XI, params_i.M, params_i.D, K_r);
+    #endif
+
+    // -- Check results --
+
+    // Simple general tests
+    for (UINT i=0; i<params_i.K; i++){
+        COMPLEX bsr = bound_states_r[i];
+        COMPLEX ncr = normconsts_r[i];
+
+        UINT is_bsr_pure_imaginary = FABS(CREAL(bsr)) < 1e-8;
+        UINT is_bsr_positive_imaginary = (CIMAG(bsr) > 0) && is_bsr_pure_imaginary;
+
+        UINT is_ncr_real = FABS(CIMAG(ncr)) < 1e-8;
+
+        if (is_bsr_pure_imaginary &&
+            is_bsr_positive_imaginary &&
+            is_ncr_real) 
+        {
+            ret_code = SUCCESS;
+        } 
+        else {
+            ret_code = FNFT_EC_TEST_FAILED;
+            break;
+        }
+    }
+    CHECK_RETCODE(ret_code, leave_fun);
+
+    // Check matching of the computed  bound states
+    // bound_states_r is sorted in ascending order, K_r is the number of found bound states by fnft_kdvv
+    // take only the last params_i.K bound states, because they are the bigger ones and more likely the ones
+    // which corresponds to the initial given bound states
+    COMPLEX * const candidate_bound_states_r_ptr = &bound_states_r[K_r-params_i.K];
+    REAL hausdorff_dist_bound_states = misc_hausdorff_dist_normed(  params_i.K, params_i.bound_states, params_i.K, 
+                                                                    candidate_bound_states_r_ptr);
+    UINT is_bsr_in_tolerance = hausdorff_dist_bound_states < err_bnd_bound_states;                                                                
+
+    COMPLEX * const candidate_normconsts_r_ptr = &normconsts_r[K_r-params_i.K];
+    REAL hausdorff_dist_normconsts = misc_hausdorff_dist_normed(params_i.K, params_i.normconsts, params_i.K, 
+                                                                candidate_normconsts_r_ptr);
+    UINT is_ncr_in_tolerance = hausdorff_dist_normconsts < err_bnd_normconst;                                                          
+
+    #ifdef DEBUG
+        printf("Number of bound_states:\n  K_r = %u\n", (unsigned int)K_r);
+        printf("Hausdorff dist bound states:\n  dist = %f\n", hausdorff_dist_bound_states);
+        printf("Hausdorff dist normconsts:\n  dist = %f\n", hausdorff_dist_normconsts);
+    #endif
+
+    UINT is_contspec_small = 1;                                                                     
+    for (UINT i=0; i<params_i.M; i++){
+        if (CABS(contspec_r[i]) > err_bnd_contspec){
+            is_contspec_small = 0;
+        }
+    }
+
+    UINT are_spurious_bound_states_small = 1;
+    for (UINT i=0; i<(K_r-params_i.K); i++){
+        if (CABS(bound_states_r[i]) > err_bnd_spurious_bound_states){
+            are_spurious_bound_states_small = 0;
+        }
+    }
+
+    if (is_bsr_in_tolerance &&
+        is_ncr_in_tolerance &&
+        is_contspec_small &&
+        are_spurious_bound_states_small) 
+    {
+        ret_code = SUCCESS;
+    } 
+    else {
+        ret_code = FNFT_EC_TEST_FAILED;
+    }
+    
+
+leave_fun:
+    free(q);
+    free(contspec_r);
+    free(bound_states_r);
+    free(normconsts_r);
+
+    if (ret_code != SUCCESS)
+        return EXIT_FAILURE;
+    else
+	    return EXIT_SUCCESS;
+}
\ No newline at end of file
diff --git a/src/private/fnft__misc.c b/src/private/fnft__misc.c
index 1dc90b6..d2b31ad 100644
--- a/src/private/fnft__misc.c
+++ b/src/private/fnft__misc.c
@@ -17,6 +17,7 @@
 * Sander Wahls (TU Delft) 2017.
 * Peter J Prins (TU Delft) 2018-2020.
 * Shrinivas Chimmalgi (TU Delft) 2019-2020.
+* Fabian Fischer (Hiwi KIT) 2026
 */
 #define FNFT_ENABLE_SHORT_NAMES
 
@@ -28,27 +29,31 @@
 void misc_print_buf(const UINT len, COMPLEX const * const buf,
                     char const * const varname)
 {
+    fnft_printf_ptr_t printf_ptr = fnft_errwarn_getprintf();
+
     UINT i;
-    printf("%s = [", varname);
+    printf_ptr("%s = [", varname);
     for (i = 0; i < len; i++) {
-        printf("%1.12e+%1.12ej", CREAL(buf[i]), CIMAG(buf[i]));
+        printf_ptr("%1.12e+%1.12ej", CREAL(buf[i]), CIMAG(buf[i]));
         if (i != len-1)
-            printf(", ");
+            printf_ptr(", ");
     }
-    printf("];\n");
+    printf_ptr("];\n");
 }
 
 void misc_print_buf_real(const UINT len, REAL const * const buf,
                     char const * const varname)
 {
+    fnft_printf_ptr_t printf_ptr = fnft_errwarn_getprintf();
+    
     UINT i;
-    printf("%s = [", varname);
+    printf_ptr("%s = [", varname);
     for (i = 0; i < len; i++) {
-        printf("%1.12e", buf[i]);
+        printf_ptr("%1.12e", buf[i]);
         if (i != len-1)
-            printf(", ");
+            printf_ptr(", ");
     }
-    printf("];\n");
+    printf_ptr("];\n");
 }
 
 REAL misc_rel_err(const UINT len, COMPLEX const * const vec_numer,
@@ -107,6 +112,38 @@ REAL misc_hausdorff_dist(const UINT lenA,
     return max_dist;
 }
 
+REAL misc_hausdorff_dist_normed(const UINT lenA,
+    COMPLEX const * const vecA, const UINT lenB,
+    COMPLEX const * const vecB)
+{
+    UINT i, j;
+    double tmp, dist, max_dist = -1.0;
+
+    for (i=0; i<lenA; i++) {
+        dist = INFINITY;
+        for (j=0; j<lenB; j++) {
+            tmp = CABS(vecA[i] - vecB[j])/CABS(vecA[i]);
+            if (tmp < dist)
+                dist = tmp;
+        }
+        if (dist > max_dist)
+            max_dist = dist;
+    }
+
+    for (j=0; j<lenB; j++) {
+        dist = INFINITY;
+        for (i=0; i<lenA; i++) {
+            tmp = CABS(vecA[i] - vecB[j])/CABS(vecB[j]);
+            if (tmp < dist)
+                dist = tmp;
+        }
+        if (dist > max_dist)
+            max_dist = dist;
+    }
+
+    return max_dist;
+}
+
 COMPLEX misc_sech(COMPLEX Z)
 {
     return 2.0 / (CEXP(Z) + CEXP(-Z));
diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c
new file mode 100644
index 0000000..b8ae38a
--- /dev/null
+++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_analytic_test.c
@@ -0,0 +1,187 @@
+/*
+* This file is part of FNFT.  
+*                                                                  
+* FNFT is free software; you can redistribute it and/or
+* modify it under the terms of the version 2 of the GNU General
+* Public License as published by the Free Software Foundation.
+*
+* FNFT is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*                                                                      
+* You should have received a copy of the GNU General Public License
+* along with this program. If not, see <http://www.gnu.org/licenses/>.
+*
+* Contributors:
+* Fabian Fischer (Hiwi KIT) 2026
+*/
+
+#define FNFT_ENABLE_SHORT_NAMES
+
+#include <stdio.h>
+
+#include "fnft__kdvv_inverse_testcases.h"
+#include "fnft__errwarn.h"
+#include "fnft__misc.h"
+
+
+/* This is a testcase with analytic reference
+ *
+ * This test case is based on the example discussed on p. 74-78 in [1]. Similar to this example
+ * is the example discussed in section 5.2 (p. 13) in [2]. Note that we use different notation and 
+ * the definition of the norming constants in [2]. Thus we have to do the following adaptions of the
+ * formulas given in [1]:
+ *
+ *   - substitute x with t and t with x
+ *   - the analytic signal needs to be multiplied by -1
+ *   - norming constants:
+ *       b(2i, x) = exp(8i*(2i)^3*x) --> b(2i, 0) = 1
+ *       b(1i, x) = -exp(8i*(1i)^3*x) --> b(1i, 0) = -1
+ *    
+ * [1] P. G. Drazin, R. S. Johnson (1989). Solitons - an introduction. Cambridge University Press
+ * [2] Prins, P. J., & Wahls, S. (2021). An accurate O(N^2) floating point algorithm for the Crum 
+ * transform of the KdV equation. Communications in Nonlinear Science and Numerical Simulation,
+ * 102, Article 105782. https://doi.org/10.1016/j.cnsns.2021.105782
+ */
+
+
+#define K 2
+
+#define QUADRATIC_ERROR_SUM_TOLERANCE 1e-27
+#define MAX_QUADRATIC_ERROR 1e-28
+
+static REAL analytic_signal(REAL t, REAL x){
+    return 12.0*(3.0 + 4.0 * COSH(2.0*t-8.0*x) + COSH(4.0*t-64.0*x))/POW((3.0*COSH(t-28.0*x) + COSH(3.0*t-36.0*x)), 2);
+}
+
+
+INT main()
+{
+    INT ret_code = SUCCESS;
+
+    COMPLEX * t_grid = NULL;
+    COMPLEX * q_1 = NULL;
+    COMPLEX * q_2 = NULL;
+    COMPLEX * analytic_q_1 = NULL;
+    COMPLEX * analytic_q_2 = NULL;
+
+    UINT D = 256;
+    REAL const T[2] = {-8.0, 12.0};
+
+    t_grid = malloc(D * sizeof(COMPLEX));
+    CHECK_NOMEM(t_grid, ret_code, leave_fun);
+    
+    COMPLEX const eps_t = (T[1] - T[0])/(D - 1);
+    
+    for (UINT n=0; n<D; n++) {
+        t_grid[n]= T[0] + n*eps_t;
+    }
+
+    // Calculating norming constants for both x-values.
+    // Attention: both norming constants arrays need to be calculated here! If there are further
+    // commands between the calculation of both norming constants arrays and the compiler optimize 
+    // drastically, there may be errors. Is is speculated, that in that case the processor 
+    // calculates the 2. norming constants with rounded intermediate results. 
+    
+    COMPLEX const bound_states[K] = { I*2.0, I*1.0 };
+    
+    REAL const x_1 = -0.1;
+    COMPLEX normconsts_1[K] = { CREAL(CEXP(8.0*I*CPOW(bound_states[0], 3)*x_1)), 
+                                -CREAL(CEXP(8.0*I*CPOW(bound_states[1], 3)*x_1)) };
+
+    REAL const x_2 = 0.5;                                
+    COMPLEX normconsts_2[K] = { CREAL(CEXP(8.0*I*CPOW(bound_states[0], 3)*x_2)), 
+                                -CREAL(CEXP(8.0*I*CPOW(bound_states[1], 3)*x_2)) };
+
+
+    // Value #1 for x                                
+    q_1 = malloc(D * sizeof(COMPLEX));
+    CHECK_NOMEM(q_1, ret_code, leave_fun);
+    
+    ret_code = fnft_kdvv_inverse(0, NULL, NULL, K, bound_states, normconsts_1, D, q_1, T, NULL);
+    CHECK_RETCODE(ret_code, leave_fun);
+    
+    analytic_q_1 = malloc(D * sizeof(COMPLEX));
+    CHECK_NOMEM(analytic_q_1, ret_code, leave_fun);
+    
+    for (UINT i=0; i < D; i++) {
+        analytic_q_1[i] = analytic_signal(t_grid[i], x_1);
+    }
+    
+    REAL error_1 = 0.0;
+    REAL max_quadratic_error_1 = 0.0;
+    REAL quadratic_error_sum_1 = 0.0;
+        
+    for (UINT i=0; i < D; i++) {
+        error_1 = CABS(CPOW((analytic_q_1[i] - q_1[i]), 2));
+        quadratic_error_sum_1 += error_1;
+        if (error_1 > max_quadratic_error_1) {max_quadratic_error_1 = error_1;}
+    }
+
+    #ifdef DEBUG
+        // misc_print_buf(K, bound_states, "bs1");
+        // misc_print_buf(D, q_1, "q_1");
+        printf("resulting max quadratic error: %e \n", max_quadratic_error_1);
+        printf("resulting quadratic error sum: %e \n", quadratic_error_sum_1);
+    #endif
+
+
+    // Value #2 for x
+    q_2 = malloc(D * sizeof(COMPLEX));
+    CHECK_NOMEM(q_2, ret_code, leave_fun);
+
+    analytic_q_2 = malloc(D * sizeof(COMPLEX));
+    CHECK_NOMEM(analytic_q_2, ret_code, leave_fun);
+    
+    ret_code = fnft_kdvv_inverse(0, NULL, NULL, K, bound_states, normconsts_2, D, q_2, T, NULL);
+    CHECK_RETCODE(ret_code, leave_fun);
+    
+    for (UINT i=0; i < D; i++) {
+        analytic_q_2[i] = analytic_signal(t_grid[i], x_2);
+    }
+
+    REAL error_2 = 0.0;
+    REAL max_quadratic_error_2 = 0.0;
+    REAL quadratic_error_sum_2 = 0.0;
+        
+    for (UINT i=0; i < D; i++) {
+        error_2 = CABS(CPOW((analytic_q_2[i] - q_2[i]), 2));
+        quadratic_error_sum_2 += error_2;
+        if (error_2 > max_quadratic_error_2) {max_quadratic_error_2 = error_2;}
+    }
+
+    #ifdef DEBUG
+        // misc_print_buf(D, q_2, "q_2");
+        printf("resulting max quadratic error: %e \n", max_quadratic_error_2);
+        printf("resulting quadratic error sum: %e \n", quadratic_error_sum_2);
+    #endif
+
+    UINT is_quadratic_error_sum_in_tolerance =  (quadratic_error_sum_1 < QUADRATIC_ERROR_SUM_TOLERANCE) &&
+                                                (quadratic_error_sum_2 < QUADRATIC_ERROR_SUM_TOLERANCE);
+    UINT is_max_quadratic_error_in_tolerance =  (max_quadratic_error_1 < MAX_QUADRATIC_ERROR) &&
+                                                (max_quadratic_error_2 < MAX_QUADRATIC_ERROR);
+
+    
+    if (is_max_quadratic_error_in_tolerance &&
+        is_quadratic_error_sum_in_tolerance) 
+    {
+        ret_code = SUCCESS;
+    } 
+    else {
+        ret_code = FNFT_EC_TEST_FAILED;
+    }
+       
+
+leave_fun:
+    free(t_grid);
+    free(q_1);
+    free(q_2);
+    free(analytic_q_1);
+    free(analytic_q_2);
+
+    if (ret_code != SUCCESS)
+        return EXIT_FAILURE;
+    else
+	    return EXIT_SUCCESS;
+}
diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c
new file mode 100644
index 0000000..bc831b1
--- /dev/null
+++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_1.c
@@ -0,0 +1,101 @@
+/*
+* This file is part of FNFT.  
+*                                                                  
+* FNFT is free software; you can redistribute it and/or
+* modify it under the terms of the version 2 of the GNU General
+* Public License as published by the Free Software Foundation.
+*
+* FNFT is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*                                                                      
+* You should have received a copy of the GNU General Public License
+* along with this program. If not, see <http://www.gnu.org/licenses/>.
+*
+* Contributors:
+* Fabian Fischer (Hiwi KIT) 2026
+*/
+
+#define FNFT_ENABLE_SHORT_NAMES
+
+#include <stdio.h>
+
+#include "fnft__kdvv_inverse_testcases.h"
+#include "fnft__errwarn.h"
+#include "fnft__misc.h"
+
+/* This is a simple test case */
+
+INT main()
+{
+    INT ret_code = SUCCESS;
+
+    COMPLEX bound_states_i[5] = {I*SQRT(5.0/2.0), I*SQRT(4.0/2.0), I*SQRT(3.0/2.0), I*SQRT(2.0/2.0), I*SQRT(1.0/2.0)};
+    COMPLEX normconsts_i[5] = {1*1e7, -1*1e-5, 1*1, -1*0.1, 1*10};
+
+    fnft_kdvv_params kdvv_parameters = {
+        .D = 256,
+        .T = {-20.0, 20.0},
+        .K = 5,
+        .bound_states = bound_states_i,
+        .normconsts = normconsts_i,
+        .M = 10,
+        .XI = {-2.0, 2.0},
+        .contspec = NULL,
+    };
+
+    kdvv_parameters.contspec = malloc(kdvv_parameters.M*sizeof(COMPLEX));
+    CHECK_NOMEM(kdvv_parameters.contspec, ret_code, leave_fun);
+    
+    REAL err_bnd_bound_states = 1.4e-3;
+    REAL err_bnd_spurious_bound_states = 1e-2;
+    REAL err_bnd_normconst = 4.5e-2;
+    REAL err_bnd_contspec = 1e-1;
+
+    
+    ret_code = kdvv_testcases_get_spectrum_of_inverse(  kdvv_parameters, err_bnd_bound_states, 
+                                                        err_bnd_spurious_bound_states, 
+                                                        err_bnd_normconst, err_bnd_contspec);
+
+    CHECK_RETCODE(ret_code, leave_fun);
+
+    
+    // Check quadratic convergence
+    kdvv_parameters.D *= 2;
+    err_bnd_bound_states /= 4;
+    err_bnd_spurious_bound_states /= 4;
+    err_bnd_normconst /= 4;
+    err_bnd_contspec /= 4;
+
+    ret_code = kdvv_testcases_get_spectrum_of_inverse(  kdvv_parameters, err_bnd_bound_states, 
+                                                        err_bnd_spurious_bound_states, 
+                                                        err_bnd_normconst, err_bnd_contspec);
+
+    CHECK_RETCODE(ret_code, leave_fun);
+
+    kdvv_parameters.D *= 2;
+    err_bnd_bound_states /= 4;
+    err_bnd_spurious_bound_states /= 4;
+    err_bnd_normconst /= 4;
+    err_bnd_contspec /= 4;
+
+    ret_code = kdvv_testcases_get_spectrum_of_inverse(  kdvv_parameters, err_bnd_bound_states, 
+                                                        err_bnd_spurious_bound_states, 
+                                                        err_bnd_normconst, err_bnd_contspec);
+
+    CHECK_RETCODE(ret_code, leave_fun);
+    
+
+leave_fun:
+    free(kdvv_parameters.contspec);    
+
+    if (ret_code != SUCCESS)
+        return EXIT_FAILURE;
+    else
+	    return EXIT_SUCCESS;
+}
+
+
+
+
diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_2.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_2.c
new file mode 100644
index 0000000..ef22dfe
--- /dev/null
+++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_2.c
@@ -0,0 +1,111 @@
+/*
+* This file is part of FNFT.  
+*                                                                  
+* FNFT is free software; you can redistribute it and/or
+* modify it under the terms of the version 2 of the GNU General
+* Public License as published by the Free Software Foundation.
+*
+* FNFT is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*                                                                      
+* You should have received a copy of the GNU General Public License
+* along with this program. If not, see <http://www.gnu.org/licenses/>.
+*
+* Contributors:
+* Fabian Fischer (Hiwi KIT) 2026
+*/
+
+#define FNFT_ENABLE_SHORT_NAMES
+
+#include <stdio.h>
+
+#include "fnft__kdvv_inverse_testcases.h"
+#include "fnft__errwarn.h"
+#include "fnft__misc.h"
+
+/* This is a testcase with many eigenvalues */
+
+#define K_I 19
+
+INT main()
+{
+    INT ret_code = SUCCESS;
+
+    COMPLEX bound_states_i[K_I] = { 40.0, 39.0, 38.0, 37.0, 36.0,
+                                    30.0, 29.0, 28.0, 27.0, 26.0,
+                                    20.0, 19.0, 18.0, 17.0, 16.0,
+                                    10.0, 9.0, 8.0, 7.0};
+
+    COMPLEX normconsts_i[K_I] = {   1e20, -1e-7, 1e5, -1e3, 1e1,
+                                    -1e0, 1e2, -1e4, 1e-6, -1e8,
+                                    1e2, -1e4, 1e6, -1e8, 1e-10,
+                                    -1e7, 1e-6, -1e5, 1e-9};
+
+    for (UINT i = 0; i<K_I; i++){
+        bound_states_i[i] = I*SQRT(bound_states_i[i]/2.0);
+    }
+
+
+    fnft_kdvv_params kdvv_parameters = {
+        .D = 512,
+        .T = {-15.0, 15.0},
+        .K = K_I,
+        .bound_states = bound_states_i,
+        .normconsts = normconsts_i,
+        .M = 10,
+        .XI = {-2.0, 2.0},
+        .contspec = NULL,
+    };
+
+    kdvv_parameters.contspec = malloc(kdvv_parameters.M*sizeof(COMPLEX));
+    CHECK_NOMEM(kdvv_parameters.contspec, ret_code, leave_fun);
+    
+    REAL err_bnd_bound_states = 2e-3;
+    REAL err_bnd_spurious_bound_states = 0.1;
+    REAL err_bnd_normconst = 1.6e-1;
+    REAL err_bnd_contspec = 2.5e-1;
+
+    
+    ret_code = kdvv_testcases_get_spectrum_of_inverse(  kdvv_parameters, err_bnd_bound_states, 
+                                                        err_bnd_spurious_bound_states, 
+                                                        err_bnd_normconst, err_bnd_contspec);
+
+    CHECK_RETCODE(ret_code, leave_fun);
+
+    
+    // Check quadratic convergence
+    kdvv_parameters.D *= 2;
+    err_bnd_bound_states /= 4;
+    err_bnd_spurious_bound_states /= 4;
+    err_bnd_normconst /= 4;
+    err_bnd_contspec /= 4;
+
+    ret_code = kdvv_testcases_get_spectrum_of_inverse(  kdvv_parameters, err_bnd_bound_states, 
+                                                        err_bnd_spurious_bound_states, 
+                                                        err_bnd_normconst, err_bnd_contspec);
+
+    CHECK_RETCODE(ret_code, leave_fun);
+
+    kdvv_parameters.D *= 2;
+    err_bnd_bound_states /= 4;
+    err_bnd_spurious_bound_states /= 4;
+    err_bnd_normconst /= 4;
+    err_bnd_contspec /= 4;
+
+    ret_code = kdvv_testcases_get_spectrum_of_inverse(  kdvv_parameters, err_bnd_bound_states, 
+                                                        err_bnd_spurious_bound_states, 
+                                                        err_bnd_normconst, err_bnd_contspec);
+
+    CHECK_RETCODE(ret_code, leave_fun);
+    
+
+leave_fun:
+    free(kdvv_parameters.contspec); 
+
+    if (ret_code != SUCCESS)
+        return EXIT_FAILURE;
+    else
+	    return EXIT_SUCCESS;
+}
diff --git a/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c
new file mode 100644
index 0000000..b825ff0
--- /dev/null
+++ b/test/fnft_kdvv_inverse/fnft_kdvv_inverse_test_3.c
@@ -0,0 +1,132 @@
+/*
+* This file is part of FNFT.  
+*                                                                  
+* FNFT is free software; you can redistribute it and/or
+* modify it under the terms of the version 2 of the GNU General
+* Public License as published by the Free Software Foundation.
+*
+* FNFT is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*                                                                      
+* You should have received a copy of the GNU General Public License
+* along with this program. If not, see <http://www.gnu.org/licenses/>.
+*
+* Contributors:
+* Fabian Fischer (Hiwi KIT) 2026
+*/
+
+#define FNFT_ENABLE_SHORT_NAMES
+
+#include <stdio.h>
+
+#include "fnft__kdvv_inverse_testcases.h"
+#include "fnft__errwarn.h"
+#include "fnft__misc.h"
+
+/* This is a testcase with an assymetric window */
+
+#define K_I 8
+
+INT main()
+{
+    INT ret_code = SUCCESS;
+
+    COMPLEX bound_states_i[K_I] = { 40.0, 35.0, 28.0, 23.0, 19.0,
+                                    16.0, 12.0, 4.0};
+
+    COMPLEX normconsts_i[K_I] = {   1e-3, -1e-17, 1e-15, -1e-20, 1e-1,
+                                    -1e-5, 1e-12, -1e-14};
+
+    for (UINT i = 0; i<K_I; i++){
+        bound_states_i[i] = I*SQRT(bound_states_i[i]/2.0);
+    }
+
+
+    fnft_kdvv_params kdvv_parameters = {
+        .D = 256,
+        .T = {-18.0, 5.0},
+        .K = K_I,
+        .bound_states = bound_states_i,
+        .normconsts = normconsts_i,
+        .M = 10,
+        .XI = {-2.0, 2.0},
+        .contspec = NULL,
+    };
+
+    kdvv_parameters.contspec = malloc(kdvv_parameters.M*sizeof(COMPLEX));
+    CHECK_NOMEM(kdvv_parameters.contspec, ret_code, leave_fun);
+
+    #ifdef DEBUG
+        printf("\n Initial parameters: \n");
+        kdvv_print_spectrum(kdvv_parameters.bound_states, kdvv_parameters.normconsts, 
+                            kdvv_parameters.contspec, kdvv_parameters.XI, 
+                            kdvv_parameters.M, kdvv_parameters.D, 
+                            kdvv_parameters.K);
+    #endif
+
+
+
+    #ifdef DEBUG
+        printf("\n -- 1. Iteration: -- \n");
+    #endif
+    
+    REAL err_bnd_bound_states = 3.8e-3;
+    REAL err_bnd_spurious_bound_states = 0.12;
+    REAL err_bnd_normconst = 2e-1;
+    REAL err_bnd_contspec = 2.5e-1;
+   
+    ret_code = kdvv_testcases_get_spectrum_of_inverse(  kdvv_parameters, err_bnd_bound_states, 
+                                                        err_bnd_spurious_bound_states, 
+                                                        err_bnd_normconst, err_bnd_contspec);
+
+    CHECK_RETCODE(ret_code, leave_fun);
+
+    
+
+    #ifdef DEBUG
+        printf("\n -- 2. Iteration: -- \n");
+    #endif
+
+    // Check quadratic convergence
+    kdvv_parameters.D *= 2;
+    err_bnd_bound_states /= 4;
+    err_bnd_spurious_bound_states /= 4;
+    err_bnd_normconst /= 4;
+    err_bnd_contspec /= 4;
+
+    ret_code = kdvv_testcases_get_spectrum_of_inverse(  kdvv_parameters, err_bnd_bound_states, 
+                                                        err_bnd_spurious_bound_states, 
+                                                        err_bnd_normconst, err_bnd_contspec);
+
+    CHECK_RETCODE(ret_code, leave_fun);
+
+
+
+    #ifdef DEBUG
+        printf("\n -- 3. Iteration: -- \n");
+    #endif
+    
+    kdvv_parameters.D *= 2;
+    err_bnd_bound_states /= 4;
+    err_bnd_spurious_bound_states /= 4;
+    err_bnd_normconst /= 4;
+    err_bnd_contspec /= 4;
+
+    ret_code = kdvv_testcases_get_spectrum_of_inverse(  kdvv_parameters, err_bnd_bound_states, 
+                                                        err_bnd_spurious_bound_states, 
+                                                        err_bnd_normconst, err_bnd_contspec);
+
+    CHECK_RETCODE(ret_code, leave_fun);
+    
+
+leave_fun:
+    free(kdvv_parameters.contspec); 
+    
+
+    if (ret_code != SUCCESS)
+        return EXIT_FAILURE;
+    else
+	    return EXIT_SUCCESS;
+}