Variable stars play a crucial role as standard candles and provide valuable insights into stellar physics. They can be modeled either through fully fledged hydrodynamical simulations or analytically as systems of coupled differential equations describing the evolution of relevant physical quantities. Typically, such equations are arrived at by simplified physical assumptions concerning the conservation laws governing stellar interiors. Here we apply a data-driven technique—sparse identification of nonlinear dynamics (SINDy)—to automatically learn governing equations from observed light curves. We apply SINDy to 3100 light curves of three different variable types from the Catalina Sky Survey. The success rate depends systematically on variable type, with possible implications for variable star classification; however, it does not obviously depend on amplitude or period. Successful models can be reduced to the generalized Lienard equation
x
̈
+
(
a
+
bx
+
c
x
̇
)
x
̇
+
x
=
0
. Members of the Lienard class of ordinary differential equations, such as the well-studied van der Pol oscillator, already saw some application to variable star modeling. For a, b = 0 the equation can be solved exactly, and it admits both periodic and nonperiodic solutions. We find a condition on the coefficients of the general equation for the presence of a limit cycle, which is also observed numerically in several instances.