Discretize PDEs that are trivially ODEs to ODEs

[hersle]

Suppose that I want to solve the heat equation $\dot{u}(t, \boldsymbol{x}) = \nabla^2 u(t, \boldsymbol{x})$ in some-dimensional infinite space. I only care about the evolution of each of its modes in Fourier space, and I have spherical symmetry (or am in 1D). In other words I want to solve this PDE:

$$\dot{u}(t, k) = -k^2 u(t, k).$$

This PDE in $(t,k)$ is trivially one independent ODE in $t$ for each $k$. Going by the (WIP) PDESystem docs I would start with something like:

using ModelingToolkit

ivs = @independent_variables t k
vars = @variables u(..)
Dt = Differential(t)
eqs = [Dt(u(t, k)) ~ -k^2 * u(t, k)]
ics = [u(0, k) ~ exp(-k^2)] # equivalent to a Gaussian in position space
domains = [t ∈ (0, Inf), k ∈ (-Inf, +Inf)]
@named heat = PDESystem(eqs, ics, domains, ivs, vars)

Then I want to discretize it to ODE(s). For example to 1 ODE where $k$ is a parameter. Or ideally to an ensemble of ODEs that can be parallellized over independent $k$ (but preferably without hardcoding values or number of $k$). How should this work?

@ChrisRackauckas and discussion in #3737, #3738 and #3745 for context.