Seems like I need to register my custom boundary function like:
However, the docs were slightly outdated so I created a PR to fix it
# ╔═╡ 1ac45f27-dc97-4a29-bb42-ceff7ab6cc90
using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets
# ╔═╡ 10dee07a-9f74-4fb5-bc2b-f700c98047f9
using Plots
# ╔═╡ c50abfc0-fb9f-11ef-31e1-3395be6ca662
begin
# Method of Manufactured Solutions: exact solution
# Parameters, variables, and derivatives
@parameters t x
@variables u(..)
Dt = Differential(t)
Dxx = Differential(x)^2
initial(x) = x == 0
@register_symbolic initial(x)
# 1D PDE and boundary conditions
eq = Dt(u(t, x)) ~ Dxx(u(t, x))
bcs = [u(0, x) ~ initial(x),
u(t, 1) ~ 0,
u(t, -1) ~ 0]
# Space and time domains
domains = [t ∈ Interval(0.0, 0.1),
x ∈ Interval(-1, 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
sol = solve(prob, Tsit5(), saveat=0.02)
# Plot results and compare with exact solution
discrete_x = sol[x]
discrete_t = sol[t]
solu = sol[u(t, x)]
plt = plot()
for i in eachindex(discrete_t)
plot!(discrete_x, solu[i, :], label="Numerical, t=$(discrete_t[i])")
end
plt
end
