public double[][] getWaveFunction(double energy) {

        // FIXME magic number
        long X_POINTS_DENSITY = 40 * 1000_000_000L; // 40 dots per 1 nm

        // Number of dots in the whole MatterStack
        int numberOfDots = (int) (getWidth() * X_POINTS_DENSITY);

        // TODO try-catch or think on long-or-int problem
        double[] waveFunction = new double[numberOfDots];
        double[] xAxis = new double[numberOfDots];
        double[] xLocal;

        int counter = 0;

        Complex[] coef; // {A, B}
        Complex pi = new Complex(0, 1);
        Complex mi = new Complex(0, -1);

        // Going through all intervals to build the wave function
        for (int k = 0; k < intervals.size(); k++) {

            MatterInterval inl = intervals.get(k);

            xLocal = inl.getGrid(X_POINTS_DENSITY);

            for (int i = 0; i < xLocal.length; i++) {

                coef = getWavefunctionCoeffs(k, energy);

                // WF = | A * exp(i*k*x) + B * exp(-i*k*x) | -- doubles not complex!
                waveFunction[i + counter] = (coef[0].times(Complex.exp(pi.times(inl.getWaveNumber(energy).times(
                        xLocal[i] - inl.getXStart()))))).add(
                        coef[1].times(Complex.exp(mi.times(inl.getWaveNumber(energy).times(xLocal[i] - inl.getXStart()))))).abs();

                xAxis[i + counter] = xLocal[i];

            }
            counter += xLocal.length;
        }
        return new double[][]{xAxis, waveFunction};
    }