FiPy set of equations diverging. Suggestions?

Question

I have tried to simplify my set of PDEs:

In this new set, p and m are my variables, where all the other values are constant. The m and p values in fc²m|m|/2DA²p are treated ass mean values on each cell. Furthermore I want to fix the inlet mass flow (m(0, t)) and the outlet pressure (p(L,t)), so I have defined it as constrains.

I have two questions: my code seems to diverge, for each time step it either gives me a very high value of m in the pipeline, while the next one gives me a very negative value and so on. Also, it seems that simply define inlet mass flow and outlet pressure gives an unreal solution.

Any thoughts and suggestions? Thank you in advance, dear community, once again.

code below:

import fipy as fi
import matplotlib.pyplot as plt
import numpy as np
import scipy as sci
from matplotlib import animation
from fipy.variables.cellVariable import CellVariable

#1. Domain
L = 101
dx = .1
mesh = fi.Grid1D(nx = L, dx=dx)

#2. Parameters values (Arbitrary) 
Lambda = 0.0001 # Friction factor
D = .4 # Pipe diameter
T = 350 # Gas Temperature
c = 380
A = 3.14 * (D **2) / 4

#3. Variables
## Pressure and mass flow
p = fi.CellVariable(name="pressure", 
                    hasOld=True, 
                    mesh=mesh, 
                    value=0.)
p.setValue(10000000.) # This is the initial condition

m = fi.CellVariable(name="mass flow", 
                    hasOld=True, 
                    mesh=mesh, 
                    value=0.)
m.setValue(300.) # This is the initial condition

#4. Boundary conditions
p.constrain (2000000., where = mesh.facesRight)
p.constrain (15., where = mesh.facesLeft)
m.constrain (500., where = mesh.facesLeft)
m.constrain (500., where = mesh.facesRight)

#5. PDE
eq1 = fi.TransientTerm(var=p) == fi.ConvectionTerm(coeff = [-(c**2) / A], var=m)
eq2 =  fi.TransientTerm(var=m) == fi.ConvectionTerm(coeff = [-A], var=p) - (Lambda*(c**2)*m.cellVolumeAverage*abs(m.cellVolumeAverage)/(2*D*A*p.cellVolumeAverage))

eqn = (eq1 & eq2)

timeStepDuration = .001
steps = 50

#This make a chart of the results
for step in range(steps):
    plt.plot(np.linspace(0, L, nx), m.value)
    #plt.xlim(0,1.1)
    #plt.ylim(0, 20)
    plt.show()
    
    eqn.sweep(dt=timeStepDuration)
    p.updateOld()
    m.updateOld()