While topological derivatives have proven useful in applications of topology optimization and inverse problems, their mathematically rigorous derivation remains an ongoing research topic, in particular in the context of nonlinear partial differential equation (PDE) constraints. We present a systematic yet formal approach for the numerical computation of topological derivatives of a large class of PDE-constrained topology optimization problems with respect to arbitrary inclusion shapes. Scalar and vector-valued as well as linear and nonlinear elliptic PDE constraints are considered in two and three space dimensions including a nonlinear elasticity model and nonlinear reaction–diffusion problems. The systematic procedure follows a Lagrangian approach for computing topological derivatives. For problems where the exact formula is known, the numerically computed values show good coincidence. Moreover, by inserting the computed values into the topological asymptotic expansion, we verify that the obtained values satisfy the expected behavior also for other, previously unknown problems, indicating the correctness of the procedure.