Static and frequency-dependent polarizabilities and first hyperpolarizabilities have been calculated for HF and Ne using full configuration interaction (FCI) and a hierarchy of coupled cluster models: coupled cluster singles (CCS), an approximate coupled cluster singles and doubles model (CC2), coupled cluster singles and doubles (CCSD), an approximate coupled cluster singles, doubles, and triples model (CC3), and coupled cluster singles, doubles, and triples (CCSDT). A previous study of BH concerning FCI benchmarking has been extended to include CC3 and static CCSDT values. Systematic improvements of the polarizabilities and the hyperpolarizabilities are found going from CCS to CCSD and from CCSD to CC3 or CCSDT. Little or no improvement of the polarizabilities and no improvement of the hyperpolarizabilities are seen when going from CCS to CC2. The CCSD results represent a significant improvement over CCS and CC2 but are again surpassed by the CC3 results which agree very well with the FCI values. The relative error for the static polarizability at the CC3 level is 0.11% for Ne and, respectively, 0.16% and 0.20% for alpha(xx) and alpha(zz) of HF. For beta(zzz) and beta(zxx) the errors are 0.50% and 1.7%, respectively. Only in the challenging case of BH does CCSDT improve the CC3 values. The dispersion for the polarizabilities and hyperpolarizabilities is predicted with increasing accuracy in the CCS-CC2-CCSD-CC3 sequence as expected from the increasing accuracy of the electronic excitation energies. For all molecules the effect of orbital relaxation has been investigated for the static properties. The inclusion of orbital relaxation gives results that are somewhat different from the unrelaxed results but are in general no improvement.