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.