@parameters t x
@parameters Dn, Dp
@variables u(..) v(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

eqs  = [Dt(u(t,x)) ~ Dn * Dxx(u(t,x)) + u(t,x)*v(t,x), 
        Dt(v(t,x)) ~ Dp * Dxx(v(t,x)) - u(t,x)*v(t,x)]
bcs = [u(0,x) ~ sin(pi*x/2),
       v(0,x) ~ sin(pi*x/2),
       u(t,0) ~ 0.0, Dx(u(t,1)) ~ 0.0,
       v(t,0) ~ 0.0, Dx(v(t,1)) ~ 0.0]

domains = [t ∈ Interval(0.0,1.0),
           x ∈ Interval(0.0,1.0)]

@named pdesys = PDESystem(eqs,bcs,domains,[t,x],[u(t,x),v(t,x)],[Dn=>0.5, Dp=>2])
discretization = MOLFiniteDifference([x=>0.1],t)
prob = discretize(pdesys,discretization) # This gives an ODEProblem since it's time-dependent

N = 10000
lb = [0.01,0.01]
ub = [2.0,2.0]
par = QuasiMonteCarlo.sample(N,lb,ub,LatinHypercubeSample())
sols1 = []
@time begin
for (Dnval, Dpval) in zip(par[1,:], par[2,:])
    newprob = remake(prob, p=[Dnval, Dpval])
    push!(sols1, solve(newprob, Tsit5()))
end
end
# 5.7 s

function prob_func(prob,i,repeat)
    remake(prob,p=par[:,i])
end

ensemble_prob = EnsembleProblem(prob, prob_func=prob_func)
@time begin
sols2 = solve(ensemble_prob, Tsit5(), trajectories=N)
end
# 91.0 s