Let (−A, B, C) be a linear system in continuous time t > 0 with input and output space C and state space H. The function φ (x) (t) = Ce −(t+2x)A B determines a Hankel integral operator) defines the tau function of (−A, B, C). Such tau functions arise in Tracy and Widom's theory of matrix models, where they describe the fundamental probability distributions of random matrix theory. Dyson considered such tau functions in the inverse spectral problem for Schrödinger's equation −f ′′ + uf = λf , and derived the formula for the potential u(x) = −2 d 2 dx 2 log τ (x) in the self-adjoint scattering case Commun. Math. Phys. 47 (1976), 171-183. This paper introduces a operator function R x that satisfies Lyapunov's equation dR x dx = −AR x − R x A and τ (x) = det(I + R x ), without assumptions of self-adjointness. When −A is sectorial, and B, C are Hilbert-Schmidt, there exists a non-commutative differential ring A of operators in H and a differential ring homomorphism ⌊ ⌋ : A → C[u, u ′ , . . .] such that u = −4⌊A⌋, which provides a substitute for the multiplication rules for Hankel operators considered by Pöppe, and McKean Cent. Eur. J. Math. 9 (2011), 205-243. The paper obtains conditions on (−A, B, C) for Schrödinger's equation with meromorphic u to be integrable by quadratures. Special results apply to the linear systems associated with scattering u, periodic u and elliptic u. The paper constructs a family of solutions to the Kadomtsev-Petviashivili differential equations, and proves that certain families of tau functions satisfy Fay's identities.