Можно с помощью пакета символьных вычислений (Symbolic toolbox):
c = 299792.458;
Omega_m = .3;
Omega_lam = .7;
H = 2.2e-18;
syms z
f = 5 * log(c * H^(-1) * (1+z) *...
int((...
(1+z)^2 * (1+Omega_m*z) - z*(2+z)*Omega_lam...
)^(-1/2), [0, z]) )+ 25;
fplot(f, [0 5])
