Hello,
You can change your code a bit like this
%% Quantum Control
format long
clear all
clc
%% Variables
time = -20000:.1:20000; % Time
n = length(time);
U = zeros(1,n);
for i = 1:n % Time
ft = time(i);
alpha = 0.001; % Transition parameter
beta = 0.001; % Transition parameter
ai = 0.2; % Initial population |1>
ai2 = 0.8; % Initial population |2>
af = 1; % Final population |1>
w0 = 0.02; % w0 = (E2 - E1)/h
mi = 6; % dipole operator
phii = pi/24; % initial relative phase
phif = pi/4; % final relative phase
E1 = 1.98; % eigen-energy 1
E2 = 2; % eigen-energy 2
H = [0 1;1 0]; % Interation hamiltonian
H0 = [E1 0;0 E2]; % hamiltonian
%% Control function
g = 1./(1+exp(-alpha.*ft));
f = ai*(1-g)+af*g;
p = 1./(1+exp(-beta.*ft));
h = phii*(1-p)+phif*p;
%% Electric field
E = alpha.*(af-ai).*exp(alpha.*ft).*sin(w0.*ft+h)./(mi.*(1+exp(alpha.*ft)).*sqrt((1-ai+(1-af).*exp(alpha.*ft)).*(ai+af.*exp(alpha.*ft))))+2.*beta.*(phif-phii).*exp(beta.*ft).*sqrt(f.*(1-f)).*cos(w0.*ft+h)/(mi.*(1-2.*f).*(1+exp(beta.*ft).^2));
%% Time evolution operator
[V, D]=eig(H);
Z = expm(H0);
Z1 = expm(D);
U(i) = exp(-1i*ft/2)*Z*exp(1i*mi*E*ft)*V*Z1*inv(V)*exp(-1i*ft/2)*Z
end
%% End
Now, U will be a vector.
