MatLab: Plotting Output Electric Field Spectrum

Question

Goal

This code combines two matrices, diagonalizes that matrix, and finds the right and left eigenvectors, right and left eigenvectors are normalized with respect to each other, then it is multiplied by an exponential function that is a function of the Eigen values and plots the output spectrum with respect to omega(w: frequency)

Code

 %%Main V2

%has an input electric field and outputs the transmitted

%define some global variables
n = 100; %step count
a = linspace(1,100,n); %our dummy parameter
p = 1; %p is the index value used of a
N = 10^8; %N is the density of particles needs to be changed 
c = 3*10^8; %m/s
k = (2*pi)/c;
%tau = linspace(0,1,n);
tau = 0; % time delay 

%the Frequency grid
w = linspace(0,1,n);
wp = linspace(0,1,n);

%%define a diagonal X0 s.t it is zero on the bounds
chi_0 = (1+1i)*(sin(5*pi*w).^2); 
%chi_0 get diagonalized in EigVect function 

%%define a 3D function same as above 
%X_nd non diagonal 
for i=1:n
   for j=1:n
     X_nd(i,j) = (1+1i)*sin(5*pi.*w(i).*wp(j)).^2;
   end 
end

%%Diagonalize Xij and Normalize EigVectors

[U_r, L, U_L] = EigVect(chi_0,X_nd,a(p),w,n);
%get the right/left eigent vectors and eigen values L
[U_rbar,U_Lbar] = Norm(U_r,U_L);
%normalize the right/left eigen vectors

%%Define the incoming electric field

E_0 = 1; %amplitude
w_0 = 0.5; %choose value from 0.3-0.7 doesn't matter for now
sigma = 0.1; 

%This for loop creates a vector of E 
for i = 1:n
    E(i) = E_0.*sin(pi*w(i))*exp(-((w(i)-w_0)/(sigma*sqrt(2))).^2);
end

%%Compute the transmitted spectrum E(omega,tau;a,N) at the exit of the sample


%exponential function of lambda
Exp_lambda = exp(N*1i*k*diag(L));
%time delay
tau_phase = exp(1i*tau.*w);

UXU = U_rbar*diag(Exp_lambda)*(U_Lbar');

%Output Electric field
E_out = UXU*(tau_phase.*E).';

hold on
plot(w,E, 'b');
plot(w,abs(E_out), 'r'); %for output plot absolute value 
legend('Input','Output')
title(['a = ',num2str(a(p)),', N = ',num2str(N),', tau = ',num2str(tau)])

%%Functions 
%Normalization of the eigenvectors w.r.t one another
function [A_bar,B_bar] = Norm(A,B)
    f = sum(B'*A,1);
    sqrt_f = sqrt(f);

    A_bar = A./sqrt_f;
    B_bar = B./sqrt_f;
end

%input a chi_0 and a non diagonal portion return right and left eigvectors
%and the diagonal eigen values

function [A_bar, L, B_bar] = EigVect(x,y,a,w,n)
    X_0 = diag(x);

    dw = 1/(n-1);

    for i = 1:n
        for j = 1:n
            X_ij(i,j) = (a*dw*y(i,j)+X_0(i,j))*sqrt(w(i)*w(j));
        end
    end
    %y is X_nd

    [A_bar, L, B_bar] = eig(X_ij);
    %returns the non-normalized right and left eigenvectors of X_ij 
    %and the diagonal eigenvalue matrix

end

I am looking for a general code review as well as any suggestions to bring the computation time down. Thank you!