Thank you!
Using u(something)
worked.
function foo(u)
a = u(t,x,θ)^2
return a
end
# Initial and boundary conditions
X = collect(range(-1, 1, step=0.1))
bcs = Vector{ModelingToolkit.Equation}(undef, length(U)+2)
for i in eachindex(X)
a = foo(u)
bcs[i] = u(0,X[i],θ) ~ -a*sin(X[i]*pi)
end