We study the Ising model two-point diagonal correlation function $ C(N,N)$ by presenting an exponential and form factor expansion in an integral representation which differs from the known expansion of Wu, McCoy, Tracy and Barouch. We extend this expansion, weighting, by powers of a variable $lambda$, the $j$-particle contributions, $ f^{(j)}_{N,N}$. The corresponding $ lambda$ extension of the two-point diagonal correlation function, $ C(N,N; lambda)$, is shown, for arbitrary $lambda$, to be a solution of the sigma form of the Painlev{e} VI equation introduced by Jimbo and Miwa. Linear differential equations for the form factors $ f^{(j)}_{N,N}$ are obtained and shown to have both a ``Russian doll nesting, and a decomposition of the differential operators as a direct sum of operators equivalent to symmetric powers of the differential operator of the elliptic integral $ E$. Each $ f^{(j)}_{N,N}$ is expressed polynomially in terms of the elliptic integrals $ E$ and $ K$. The scaling limit of these differential operators breaks the direct sum structure but not the ``Russian doll structure. The previous $ lambda$-extensions, $ C(N,N; lambda)$ are, for singled-out values $ lambda= cos(pi m/n)$ ($m, n$ integers), also solutions of linear differential equations. These solutions of Painleve VI are actually algebraic functions, being associated with modular curves.