Use a symbolic function (in R2012a and later versions):
syms y1 y2 b1x b2x s11 s22 s12 s k t
q = s11*(b1x - y1)^2 + s22*(b2x - y2)^2 + s12*(b1x - y1)*(b2x - y2);
fun(y1) = k*y1 * exp(-q/(2*t));
f = int(fun,y1)
f = simplify(f, 'Steps',20)
producing:
f(y1) =
- (k*t*exp(-(s11*b1x^2 + s12*b1x*b2x - 2*s11*b1x*y1 - s12*b1x*y2 + s22*b2x^2 - s12*b2x*y1 - 2*s22*b2x*y2 + s11*y1^2 + s12*y1*y2 + s22*y2^2)/(2*t)))/s11 - (2^(1/2)*k*pi^(1/2)*erf((2^(1/2)*(b1x*s11*2i + b2x*s12*1i - s11*y1*2i - s12*y2*1i))/(4*t*(-s11/t)^(1/2)))*exp(-((b2x - y2)^2*(- s12^2 + 4*s11*s22))/(8*s11*t))*(2*b1x*s11 + b2x*s12 - s12*y2)*1i)/(4*s11*(-s11/t)^(1/2))
