$(1)\ E(\hat{Y_h}) = E(b_0 + b_1X_h) = \beta_0 +\beta_1X_h$
$(2)\ var(\hat{Y_h}) = var(b_0 + b_1X_h)$
An alternate (but equivalent) version of the regression model can be written as:
$Y_i = \beta_0X_0 + \beta_1X_1 + \epsilon_i$
This model associates an X variable with each coefficient $(where X_0 = 1)$
Al alternate modification is to use the deviation $X_i -\bar{X}$ rather than $X_i$
So the model can be written as:
$Y_i = \beta_0^* + \beta_1(X_i - \bar{X}) + \epsilon_i$
where $(3)\ \beta_0^* = \beta_0 + \beta_1\bar{X}$
These models can be used interchangably.
We know from the normal equations:
$\Sigma Y_i = nb_0 + b_1\Sigma X_i$
solving for $b_0$
$(4)\ b_0 = \bar{Y} - b_1\bar{X}$
So substituting from (3) and (4):
$b_0^* = b_0 + b_1\bar{X} = (\bar{Y} - b_1\bar{X}) + b_1\bar{X} = \bar{Y}$
$(5)\ var(\hat{Y_h}) = var(b_0 + b_1X_h) = var(\bar{Y} + b_1(X_h - \bar{X}))$
using:
$var(\bar{Y}) = \frac{\sigma^2}{n}$
$var(aX) = a^2var(X)$
and
$var(X + Y) = Var(X) + Var(Y) + 2Cov(X,Y)$
So:
= $var(\bar{Y}) +(X_h - \bar{X})^2var(b_1) + 2(X_h-\bar{X})cov(\bar{Y},b_1)$
we use the fact that $Cov(\bar{Y},b_1) = 0$ due to the i.i.d assumption on $\epsilon_i$
$= \frac{\sigma^2}{n} + (X_h-\bar{X})^2\frac{\sigma^2}{\Sigma(X_i-\bar{X})^2}$
$= \sigma^2[\frac{1}{n} + \frac{(X_h - \bar{X})^2}{\Sigma(X_i - \bar{X})^2}]$