Interpolation in space not possible if solution is only saved at the end.

[Qfl3x]

If the solution is constrained to the last point in time, interpolation will no be possible.

using DifferentialEquations, ModelingToolkit, MethodOfLines, DomainSets
# Method of Manufactured Solutions: exact solution
u_exact = (x,t) -> exp.(-t) * cos.(x)

# Parameters, variables, and derivatives
@parameters t x
@variables u(..)
Dt = Differential(t)
Dxx = Differential(x)^2

# 1D PDE and boundary conditions
eq  = Dt(u(t, x)) ~ Dxx(u(t, x))
bcs = [u(0, x) ~ cos(x),
        u(t, 0) ~ exp(-t),
        u(t, 1) ~ exp(-t) * cos(1)]

# Space and time domains
domains = [t ∈ Interval(0.0, 1.0),
           x ∈ Interval(0.0, 1.0)]

# PDE system
@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)])

# Method of lines discretization
dx = 0.1
order = 2
discretization = MOLFiniteDifference([x => dx], t)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)

# Solve ODE problem
using OrdinaryDiffEq
sol = solve(prob, Tsit5(), save_on=false, save_start=false);

sol(1.0, 0.3) #Last time step; crashes

[xtalax]

This is a side effect of interpolations.jl unfortunately

[ChrisRackauckas]

Interpolations.jl isn't used?

We can probably fix this though.

[xtalax]

Yes it is, this provides the multi dimensional interp in the solution interface. I can strip the Array of the short dim, and warn.

[ChrisRackauckas]

But we know the polynomial, why use Interpolations?