  We study the critical Ising model with free boundary conditions on finite
domains in $\mathbb{Z}^d$ with $d\geq4$. Under the assumption, so far only
proved completely for high $d$, that the critical infinite volume two-point
function is of order $|x-y|^{-(d-2)}$ for large $|x-y|$, we prove the same is
valid on large finite cubes with free boundary conditions, as long as $x, y$
are not too close to the boundary. This confirms a numerical prediction in the
physics literature by showing that the critical susceptibility in a finite
domain of linear size $L$ with free boundary conditions is of order $L^2$ as
$L\rightarrow\infty$. We also prove that the scaling limit of the near-critical
(small external field) Ising magnetization field with free boundary conditions
is Gaussian with the same covariance as the critical scaling limit, and thus
the correlations do not decay exponentially. This is very different from the
situation in low $d$ or the expected behavior in high $d$ with bulk boundary
conditions.
