Error using surf, Z must be a matrix, not a scalar or vector to surf(deltaT,deltaR,abs(Err))

Question

[alpha]=MIMO_elliptical_parameters(30,90,90,10,91,180,0,'NMMEA',1)
[alpha]=MIMO_elliptical_parameters(N,BetaT,BetaR,kappa,f_max,alpha_0,alphav,method,PLOT)
if nargin<9
    PLOT=0;
end
% Converting computation parameters from degrees to radians
BetaT=2*pi*BetaT/360;
BetaR=2*pi*BetaR/360;
alpha_0=2*pi*alpha_0/360;
alphav=2*pi*alphav/360;
save data_elliptical.mat N BetaT BetaR kappa f_max alpha_0 alphav;

% Computation of initial values for the parameter alpha usng EMEDS
if all(upper(method)=='EMEDS')
    alpha=zeros(N,1);
    for n=1:N
        alpha_0_r=0.5*pi/N;
        alpha(n)=(1/N)*2*pi*(n-0.5)+alpha_0_r;
    end
    
% Computation of initial values for the parameter alpha using Modified
% Extended Method of Exact Doppler Spread (MEMEDS)
elseif all(upper(method)=='MEMED')
    alpha=zeros(N,1);
    for n=1:N
        alpha_0_r=pi/N;
        alpha(n)=((1/N)*2*pi*n)+alpha_0_r;
    end
        
% Computation of initial values for the parameter alpha using MMEA
elseif all(upper(method)=='MMEAS')
    alpha_n=linspace(-pi,pi,400*N)
    alpha=zeros(N,1);
    for n=1:N
        alpha(n)=trapz(alpha_n,Mises_PDF(alpha_n,kappa,alpha_0))-((1/N)*(n-1/4));
    end
else all(upper(method)=='NMMEA')
    alpha_n=linspace(-pi,pi,400*N);
    alpha=zeros(N,1);
    for n=1:N
        alpha(n)=trapz(alpha_n,Mises_PDF(alpha_n,kappa,alpha_0))-((1/N)*(n-0.55));
end

% testing the correctness of the simulation model by determining the time
% ACF for the simulation model as well as the reference model.
tau=linspace(0,N/(4*f_max),6*N);
alphaMS=linspace(-pi,pi,4*N);
for k=1:length(tau)
    ACF(k)=trapz(alphaMS,exp(-sqrt(-1)*2*pi*f_max*cos(alphaMS-alphav)*tau(k)).*Mises_PDF(alphaMS,kappa,alpha_0));
    ACF_teld(k)=sum(exp(-sqrt(-1).*2*pi.*f_max.*cos(alpha-alphav).*tau(k)))./length(alpha);
end

% testing the correlation of the simulation model by determining the
% 2D-space CCF for the simulation model as wel as the reference model
deltaT=linspace(0,N/5,2*N);
deltaR=linspace(0,N/5,2*N);
alphaR=linspace(-pi,pi,4*N);
CCF=zeros(length(deltaT),length(deltaR));
CCF_teld=zeros(length(deltaT),length(deltaR));
for q=1:length(deltaT)
    for g=1:length(deltaR)
        aa=exp(-sqrt(-1).*2*pi.*deltaT(q).*cos(gen_alphaT(alphaR)-BetaT));
        bb=exp(-sqrt(-1).*2*pi.*deltaR(g).*cos(alphaR-BetaR));
        CCF(q,g)=trapz(alphaR,aa.*bb.*Mises_PDF(alphaR,kappa,alpha_0));
        aa_simu=exp(-sqrt(-1).*2*pi.*deltaT(q).*cos(gen_alphaT(alphaR)-BetaT));
        bb_simu=exp(-sqrt(-1).*2*pi.*deltaR(g).*cos(alphaR-BetaR));
        CCF_teld(q,g)=(sum(aa_simu.*bb_simu)/length(alpha));
    end
deltaT_max=deltaT(length(deltaT));
deltaR_max=deltaR(length(deltaR));
if kappa==0
    Err=norm(CCF-CCF_teld,2);
else
    Err=(1/deltaT_max*deltaR_max)*trapz(deltaT,trapz(deltaR,abs(CCF-CCF_teld).^2))^1/2;

if PLOT==1
    subplot(2,2,1)
    grid on;
    plot(tau*f_max,real(ACF),'g*');
    hold on;
    plot(tau*f_max,real(ACF_teld),'r*')
    xlim([0 N/2])
    ylim([-1 1])
    set(0,'DefaultAxesFontSize',16);
    legend({['Reference model'],['Simulation model with N=',num2str(N)]},'FontName','Arial','Location','NorthEast','FontSize',16,'Interpreter','latex');
    xlabel({'$\tau$ . $f_{max}$'},'FontName','Arial','FontSize',20,'Interpreter','latex');
    ylabel({'Time ACF'},'FontName','Arial','FontSize',20,'Interpreter','latex');
    
    subplot(2,2,2)
    surf(deltaR,deltaT,(real(CCF)));
    zlim([-.5 1])
    ylim([0 5])
    xlim([0 5])
    set(gca,'YDir','reverse');
    set(0,'DefaultAxesFontSize',18);
    title({'Reference model'},'FontName','times','FontSize',20,'Interpreter','Latex');
    zlabel({'2D space CCF,$\rm\mid{{\rho_{11,22}}\mbox(\delta_T,\delta_R)}\mid$'},'FontName','times','FontSize',20,'Interpreter','latex');
    xlabel({'$\delta_{R}$ / $\lambda_0$'},'FontName','times','FontSize',20,'Rotation',12,'Interpreter','latex');
    ylabel({'$\delta_{T}$ / $\lambda_0$'},'FontName','times','FontSize',20,'Rotation',-20,'Interpreter','latex');
    
    subplot(2,2,3)
    surf(deltaR,deltaT,real(CCF_teld));
    zlim([-.5 1])
    ylim([0 5])
    xlim([0 5])
    set(gca,'YDir','reverse');
    set(0,'DefaultAxesFontSize',18);
    title({'Simulation model'},'FontName','times','FontSize',20,'Interpreter','latex');
    zlabel({'2D space CCF,$\rm{\tilde{\rho}_{11,22}\mbox(\delta_T,\delta_R)}$'},'FontName','times','FontSize',20,'Interpreter','latex');
    xlabel({'$\delta_{R}$ / $\lambda_0$'},'FontName','times','FontSize',20,'Rotation',12,'Interpreter','latex');
    ylabel({'$\delta_{T}$ / $\lambda_0$'},'FontName','times','FontSize',20,'Rotation',-20,'Interpreter','latex');
    
    subplot(2,2,4)
    surf(deltaT,deltaR,abs(Err));
    zlim([0 1])
    ylim([0 5])
    xlim([0 5])
    set(gca,'YDir','reverse');
    set(0,'DefaultAxesFontSize',18);
    title({'Simulation model'},'FontName','times','FontSize',20,'Interpreter','latex');
    zlabel({'Absolute error, $\rm{e_T\mbox(\delta_T,\delta_R)}$'},'FontName','times', 'FontSize',20,'Interpreter','latex');
    xlabel({'$\delta_{R}$ / $\lambda_0$'},'FontName','times','FontSize',20,'Rotation',12,'Interpreter','latex');
    ylabel({'$\delta_{T}$ / $\lambda_0$'},'FontName','times','FontSize',20,'Rotation',-20,'Interpreter','latex');
end
    end
end
end
end

I am trying to model an MIMO channe for wireless communication, a tpical example from Matthias Patzold book on Wireless Channel(Elliptical shape) I tried to run the Matlab codes and I could plot three (3) graphs correctly, except the fourth graph involving "surf(deltaT, deltaR,abs(Err))". That Z must be a matrix, I have tried plenty of methods to no avail. Kindly assist.