This paper presents an approach for deriving the continuum sensitivity of superconducting systems operating at critical current densities and an optimization method based on the continuum sensitivity. In the sensitivity problem, the superconducting systems is represented by a variational state equation, wherein the magnetic permeability depends on the magnetic field, which is transformed from a state equation with a field-dependent source. The design sensitivity is derived using the material derivative concept of continuum mechanics and the adjoint variable method. The adjoint system has a material property represented as a symmetric tensor that contains the sensitivity of the current density with respect to the magnetic field. The design sensitivity is represented in the analytical form of a surface integral on the interface between the superconducting material and its surroundings, which depends on the sensitivity of the current density. The optimization scheme is constructed based on the continuum design sensitivity. In the design optimization, the level set method is used to express the shape variation of the superconducting materials. The numerical example of infinite solenoids demonstrates that the design sensitivity provides an accurate design solution considering the critical current condition. In addition, the design example of a magnetic resonance imaging solenoid shows that the derived design sensitivity has the inherent ability for attaining the compact design by treating the input current of a superconducting system as a critical condition.