Solve waterhammer PDE in numpy/scipy

Question

I’ve been trying to solve the water hammer PDE’s from the Maple example linked below in python (numpy/scipy). I’m getting very unstable results. Can anyone see my mistake? Guessing something is wrong with the boundary conditions.

https://www.maplesoft.com/support/help/view.aspx?path=applications/WaterHammer

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

## Parameters
Dia = 0.1
V = 14.19058741 # Stead state
p = 1000 # Liquid density
u = 0.001 # Viscosity
L = 25 
e = 0.0001 # Roughness
Psource = 0.510E6
thick = 0.001
E= 7010*10**9
K=20010E6
Vsteady= 14.19058741
Ks = 1/((1/K)+(Dia/E*thick))

# Darcy-Weisbach
def Friction(V):
    Rey = ((Dia*V*p)/u)
    fL = 64/Rey
    fT = 1/((1.8*np.log10((6.9/Rey) + (e/(3.7*Dia))**1.11))**2)
    
    if Rey >= 0 and Rey < 2000:
        return fL
    
    if Rey >= 2000 and Rey<4000:
        return fL + ((fT-fL)*(Rey-2000))/(4000-2000)
     
    if Rey >= 4000:
        return fT
    
    return 0

def model(D, t):
    V = D[:N]
    P = D[N:]
    
    dVdt = np.zeros(N)
    for i in range(1, len(dVdt)-1):
        dVdt[i] = -(1/p)*((P[i+1]-P[i-1])/2*dx)-((Friction(np.abs(V[i]))*(np.abs(V[i])**2))/(2*Dia))
    
    dPdt = np.zeros(N)
    for i in range(1, len(dPdt)-1):
        dPdt[i] = -((V[i+1]-V[i-1])/(2*dx))*Ks
        
    if t < 2:
        dVdt[29] = 0
    else:
        dVdt[29] = -1
     
    dPdt[29] = 0
    dVdt[0] = dVdt[1]
           
    return np.append(dVdt,dPdt)


N = 30
x = np.linspace(0, L, N)
dx = x[1] - x[0]

## Initial conditions
Vi_0 = np.ones(N)*Vsteady
Pi_0 = np.arange(N)
for i in Pi_0:
    Pi_0[i] = Psource - (i*dx/L)*Psource

# initial condition
y0 = np.append(Vi_0, Pi_0)

# time points
t = np.linspace(0,3,10000)

# solve ODE
y = odeint(model,y0,t)

Vr = y[:,0:N]
Pr = y[:,N:]

plt.plot(t,Pr[:,5])