Matlab: Radix 2 FFT with 8 Channels and High Number of Inputs
Your problem is probably with twiddle factor adjustment. Below is a working code for 8 channel implementation of FFT (where the last 3 stages are explicitly computed, but not in a butterfly structure). Its derived based on the recursive structure of FFT. Where the splitting algorithm of N-point DFT into 2 N/2 point DFTs is applied for three consequtive stages.

Mathematical definition of the solution is as follows: Assume that $N$-point DFT $X[k]$ is split into two $N/2$ point DFTS according to the even and odd partitioning as usual:
$$ X[k] = X_e[k] + W_N^k X_o[k] ~~~,~~~\text{ for } k=0,1,...,N-1 $$

where the twiddle factor is $W_N = e^{-j\frac{2\pi}{N}}$ and note that the partial DFTs $X_e[k]$ and $X_o[k]$ are of length $N/2$and periodic. Then we can further decompose those $X_e[k]$ and $X_o[k]$ recursively to reach

$$\begin{align}
X_e[k] &amp;= X_{ee}[k] + W_{N/2}^k X_{eo}[k] ~~~,~~~\text{ for } k=0,1,...,N/2-1 \\
X_o[k] &amp;= X_{oe}[k] + W_{N/2}^k X_{oo}[k] ~~~,~~~\text{ for } k=0,1,...,N/2-1 \\
\end{align} $$

And further split those $N/4$ point DFTs into $N/8$ points DFTS which will be the channel signal DFTS. 

$$\begin{align}
X_{ee}[k] &amp;= X_{eee}[k] + W_{N/4}^k X_{eeo}[k] ~~~,~~~\text{ for } k=0,1,...,N/4-1 \\
X_{eo}[k] &amp;= X_{eoe}[k] + W_{N/4}^k X_{eoo}[k] ~~~,~~~\text{ for } k=0,1,...,N/4-1 \\
X_{oe}[k] &amp;= X_{oee}[k] + W_{N/4}^k X_{oeo}[k] ~~~,~~~\text{ for } k=0,1,...,N/4-1 \\
X_{oo}[k] &amp;= X_{ooe}[k] + W_{N/4}^k X_{ooo}[k] ~~~,~~~\text{ for } k=0,1,...,N/4-1 \\
\end{align} $$

Below is an implementation of the above splitting algorithm. Note that your definition of the signals X0_0,X_0_1...,X0_7 is actually correct according to this splitting mode. But the butterfly computations seems not right.

clc;clear all; clear all


% S0 - Generate N point input data sequence x[n] of a sine wave.
% --------------------------------------------------------------
N           = 128;          % Length of signal
M           = 8;            % Number of channels    
x           = sin(2*pi*0.123*[0:N-1]);

% S1 - obtain the even-odd partition signals :
% --------------------------------------------
xe = x(1:2:N);              % EVEN part of N-point x[n]
xo = x(2:2:N);              % ODD  part of N-poibt x[n]

xee  = xe(1:2:N/2);         % EVEN part of N/2-point xe[n]
xeo  = xe(2:2:N/2);         % ODD  part of N/2-point xe[n]
xoe  = xo(1:2:N/2);         % EVEN part of N/2-point xo[n]
xoo  = xo(2:2:N/2);         % ODD  part of N/2-point xo[n]


xeee = xee(1:2:N/4);         % EVEN part of N/4-point xee[n]
xeeo = xee(2:2:N/4);         % ODD  part of N/4-point xee[n]
xeoe = xeo(1:2:N/4);         % EVEN part of N/4-point xeo[n]
xeoo = xeo(2:2:N/4);         % ODD  part of N/4-point xeo[n]
xoee = xoe(1:2:N/4);         % EVEN part of N/4-point xoe[n]
xoeo = xoe(2:2:N/4);         % ODD  part of N/4-point xoe[n]
xooe = xoo(1:2:N/4);         % EVEN part of N/4-point xoo[n]
xooo = xoo(2:2:N/4);         % ODD  part of N/4-point xoo[n]

%S1 - Compute 8 of those N/8 point FFTs for the eee-ooo signals :
% ---------------------------------------------------------------
Xeee = fft(xeee,N/8);
Xeeo = fft(xeeo,N/8);
Xeoe = fft(xeoe,N/8);
Xeoo = fft(xeoo,N/8);
Xoee = fft(xoee,N/8);
Xoeo = fft(xoeo,N/8);
Xooe = fft(xooe,N/8);
Xooo = fft(xooo,N/8);

% S1 - obtain N/4 point DFT's from 8 , N/8 point FFTs of eee-ooo signals :
% ------------------------------------------------------------------------
W4 = exp(-1j*2*pi*(4/N));      % twiddle factor for the 3rd stage (Wn/4)
Xee = [Xeee, Xeee]  + (W4.^[0:N/4-1]).*[Xeeo, Xeeo];
Xeo = [Xeoe, Xeoe]  + (W4.^[0:N/4-1]).*[Xeoo, Xeoo];
Xoe = [Xoee, Xoee]  + (W4.^[0:N/4-1]).*[Xoeo, Xoeo];
Xoo = [Xooe, Xooe]  + (W4.^[0:N/4-1]).*[Xooo, Xooo];

% S2 - obtain N/2 point DFT's from N/4 point FFTs of ee-oo signals :
% ------------------------------------------------------------------
W2 = exp(-1j*2*pi*(2/N));      % twiddle factor for the 2rd stage (Wn/2)
Xe = [Xee, Xee] + (W2.^[0:N/2-1]).*[Xeo, Xeo];
Xo = [Xoe, Xoe] + (W2.^[0:N/2-1]).*[Xoo, Xoo];

% S3 - obtain N point final DFT from two N/2 point DFTs of e-o signals :
% ----------------------------------------------------------------------
Wn = exp(-1j*2*pi/N);      
X = [Xe, Xe] + (Wn.^[0:N-1]).*[Xo, Xo]; 


% SX - Display results
% --------------------
figure,stem(abs(X))
figure,stem(abs(fft(x,N)),'g');