Question regarding implementation of Multivariate Generalized Linear Model

I am trying to implement Multivariate General Linear Model in Matlab, which is implemented in R.

When I put in some random input in both R and my implementation, I get different answers and I am not sure how to proceed.

I know where I am doing it wrong, but I cannot fix it since I do not know the details of the theory.

It would be really great if anyone can tell me what I am doing wrong here.

Thanks.

if true
% MGLM using Wilk's Lambda
% Model: Y = XB + E
%
% Input
%     Y: response multivariate variable, n x m
%     X: model matrix, n x (k+1)
%     
% Output
%     stats.F: F-stats
%     stats.p: p-value
%
clc; clear;
%%Variables
mu1= [0 1 2; 0 2 2; 0 3 2;0 4 3]; % group 1
mu2=[2 0 1; 0 0 1;1 0 1;3 0 1]; % group 2
Y = [mu1;mu2];  % response variable
X = [1 1 1 1 0 0 0 0]'; % model 
%%dimensions
n=size(Y,1);
k=size(X,2); % Constant
p=size(X,2);
%%Computing sum-of-squares
temp = inv([ones(n,1) X]'*[ones(n,1) X])*[ones(n,1) X]'*Y;
intercept = temp(1,:);  % intercept
eB = temp(2,:); % B estimation
eY = X*eB;  % Y estimation
eE = Y - repmat(intercept,n,1) - eY;    % residual
Y2 = Y - repmat(intercept,n,1);
SSP_R = eE'*eE; % SSP residual matrix
y_mean = mean(Y)';  % mean of y
SSP_Reg = eY' * eY - n*y_mean*y_mean';
% Where I am doing it wrong....
L = ones(n,1);  % Hypothesis matrix L 
SSP_H = (eB'*L')*pinv(L*pinv(X'*X)*L')*(L*eB); 
%%Wilk's lambda 
WL = det(SSP_R) / det(SSP_H + SSP_R);
%%F-stat and p-value approximation 
m = size(SSP_R,1); 
s = p; 
r = n-k-1 - (m-s+1)/2;
u = (m*s-2)/4;
if m^2 + s^2 - 5 > 0
  t = sqrt((m^2 * s^2 - 4) / (m^2 + s^2 - 5));
else
  t = 0;
end
F = (1-WL^(1/t))/ WL^(1/t) * (r*t - 2*u)/(m*s);
pval = 1-fcdf(F,m*s, r*t-2*u);
%%Output
stats.F = F;
stats.p = pval;
display(stats);
end