There are different ways to do this.
This approach simply concatenates them as the third dimension of a
(4x4xnumel(B))
array —
d0=0.2;%meV delta 0
d1=0.18;%#delta 1
d2=0.05;%#delta 2
%#the electron and hole effective g factor in the z and x directions.
ghx=-0.35 ;
ghz=-2.2;
gex=-0.65;
gez=-0.8;
alpha=0.02;%#meV/T^2 diamagnetic shift coefficient,
ga=0.0 ;
gb = ga;% #斯特恩的交互gth with the exciton state with same polarization.
theta=pi/3 ;
wx=0.01;
wd=wx-d0;
wa=wx+0.5;
wb=wa;
wl=wa;
G0a=0.05;%# is the maximum value of energy of the pulse
G0b=G0a;
tc=30 ;%# the time of the center of the pulse
tal=10 ;%#is the Gaussian rms width
mub=0.0579;%# Bohr magneton
B=0:0.1:10;
V = NaN(4,4,numel(B));% Preallocate
D = V;% Preallocate
fork = 1:numel(B)
bp=mub*B(k).*sin(theta)*(gez+ghz)*0.5 ;%bp
bm=mub*B(k).*sin(theta)*(gez-ghz)*0.5;%bm
be=mub*B(k).*cos(theta)*gex*0.5;%be
bh=mub*B(k).*cos(theta)*ghx*0.5;%bh
H1=[wx,d1/2,0,0;d1/2,wx,0,0;0,0,wd,d2/2;0,0,d2/2,wd];
H2=[bp+alpha*B(k).^2,0,be,bh;0,-bp+alpha*B(k).^2,bh,be;be,bh,-bm+alpha*B(k).^2,0;bh,be,0,bm+alpha*B(k).^2];
HT=H1+H2;
[V(:,:,k),D(:,:,k)] = eig(HT);
end
Another option would be to store them in cell arrays.
.
