  In ${\cal N}=1$ SQED with $N_f$ flavors regularized by higher derivatives we
obtain the four-loop beta function using a method based on calculating vacuum
supergraphs. For this purpose we use a special C++ program which obtain
contributions to the $\beta$-function from supergraphs without external legs in
the form of integrals of double total derivatives. After that the result was
compared with the three-loop anomalous dimension calculated earlier. We
explicitly check the NSVZ relation in this order.
