Try this:
gamma_solve = @(gamma_t,gamma_0) ei(gamma_t./gamma_0)-(gamma_0/9.022);
gamma=0:2:20;
for k1 = 1:numel(gamma)
gamma0=10.^(gamma(k1)/10);
GammaT(k1) = fsolve(@(gamma_t) gamma_solve(gamma_t,gamma0), 1);
end
figure(1)
plot(gamma, GammaT)
grid
xlabel('\bf\gamma\rm')
ylabel('\bf\gamma\rm_t')
