Hi

I’m trying to recreate this Maple water hammer example in Julia. Link to Maple example below. I’m not able to get anything close to stable results. Some help here would be greatly appreciated. (I know I’ve made some simplifications, but can’t see that this should effect the results). Pretty new to Julia so there is probably a lot of improvments to be made.

```
using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets, Plots
@parameters x y t
@variables P(..) V(..)
Dt = Differential(t)
Dx = Differential(x)
friction(V) = 64/((0.1*V*1000)/0.001)
x_min = y_min = t_min = 0.0
x_max = 25.0
t_max = 3.0
α = 10.
Ks = 200000000
p = 1000
Psource = 755147.8914031894
u0(x,y,t) = 22(y*(1-y))^(3/2)
v0(x,y,t) = 27(x*(1-x))^(3/2)
P0(x,t) = ((25-x)/25)*Psource
V0(x,t) = (14.19058741*(t<=2)) +(14.19058741^(-70*(t-2)) * (t>2))
eq = [Dt(V(x,t)) + (1/p)*Dx(P(x,t)) + ((0.03*abs(V(x,t))*V(x,t))/(2*0.1)) ~ 0,
Dx(V(x,t)) + (1/Ks)+Dt(P(x,t)) ~ 0]
domains = [x ∈ Interval(x_min, x_max),
t ∈ Interval(t_min, t_max)]
# Periodic BCs
bcs = [V(x,0) ~ 14.19058741, # t=0
V(x_max,t) ~ V0(x,t),
P(x,0) ~ P0(x,t),
P(0,t) ~ Psource]
# P(x_max,t) ~ 0]
@named pdesys = PDESystem(eq,bcs,domains,[x,t],[P(x,t),V(x,t)])
N = 32
order = 2 # This may be increased to improve accuracy of some schemes
# Integers for x and y are interpreted as number of points. Use a Float to directtly specify stepsizes dx and dy.
discretization = MOLFiniteDifference([x=>N], t, approx_order=order)
# Convert the PDE problem into an ODE problem
println("Discretization:")
@time prob = discretize(pdesys,discretization)
println("Solve:")
@time sol = solve(prob, TRBDF2(), saveat=0.1)
plot(sol[V(x,t)][:,26]) ## Plotting V at time 26
```