In this paper, we present an algorithmic approach to the construction of Lyapunov functions for infinite-dimensional systems. This paper unifies and significantly extends many previous results which have appeared in conference and journal format. The unifying principle is that any linear parametrization of operators in Hilbert space can be used to construct an LMI parametrization of positive operators via squared representations. For linear systems, we get positive linear operators and hence quadratic Lyapunov functions. For nonlinear systems, we get nonlinear operators and hence non-quadratic Lyapunov functions. Special cases of these results include positive operators defined by multipliers and kernels which are: polynomial; piecewise-polynomial; or semi-separable and apply to systems with delay; multiple spatial domains; or mixed boundary conditions. We also introduce a set of efficient software tools for creating these functionals. Finally, we illustrate the approach with numerical examples.