This is not a complete answer, but I think it sheds some light on the problem.
First, let's simplify our notation the following way:
$$c=\sqrt{\frac{E}{\rho}}$$
$$\frac{1}{A}\frac{dA}{dx}=B(x)$$
Thus, we are interested in the solutions to the equation:
$$\frac{\partial^2 u}{\partial x^2}-\frac{1}{c^2} \frac{\partial^2 u}{\partial t^2}+B(x) \frac{\partial u}{\partial x}=0 \tag{1}$$
Let us consider the simple case of piecewise constant $B$, which still makes it possible to model variable cross-sections with any desired accuracy.
So, without the loss of generality we consider the solution in some region where $B$ is constant. Let us divide (1) by $B^2$ and set:
$$\chi=B x$$
$$\tau=B c t$$
We obtain a new equation:
$$\frac{\partial^2 u}{\partial \chi^2}-\frac{\partial^2 u}{\partial \tau^2}+\frac{\partial u}{\partial \chi}=0 \tag{2}$$
It's a particular case of telegraph equation and can be reduced to Klein-Gordon equation.
However, first, let's look at the form of equation and the variables. A simple wave equation has solutions in the form:
$$u_o(x,t)=f(x - vt)=f(\phi)$$
Where the phase velocity: $$v= \frac{\partial \phi}{\partial t}$$
In our case we could write (without proof so far!):
$$\phi = \chi-\tau = B x-Bct=B(x-ct)$$
Which gives us phase velocity $v=c$ which is independent of $B$, i.e. on any interval with constant $B$ the velocity would be the same, and we can take the limit to arbitrary $B(x)$, keeping the velocity constant for every $x$, just as the OP stated.
But! This argument is not rigorous. I'm not at all certain this kind of solution applies here, see below.
Now as an addition, let's try actually solving the equation (2). The link provided above recommends the following substitution:
$$u(\chi,\tau)=e^{-\chi/2} w(\chi,\tau)$$
Which reduces the equation (2) to a familiar Klein-Gordon equation:
$$\frac{\partial^2 w}{\partial \chi^2}-\frac{\partial^2 w}{\partial \tau^2}-\frac{1}{4} w=0 \tag{3}$$
If you carefully look at the section "Particular solutions" in the linked pdf, you will notice that the general form of the solutions seems familiar, however some additional conditions still apply, which makes determining the phase and the phase velocity quite difficult.
In other words, I'm not sure what the phase is supposed to be in this general case, and as a result, how to determine the phase velocity.
One suggestion I have is to solve the equation for particular $A(x)$ numerically and then see how the solution behaves.
Update
We, of course, can directly use Fourier method on the equation (1):
$$u=X(x)T(t)$$
Then the equation transforms to:
$$X''T-\frac{1}{c^2} XT''+B(x)X'T=0$$
$$\frac{X''}{X}T+B(x)\frac{X'}{X}=\frac{1}{c^2} \frac{T''}{T}$$
Since we have function of $x$ on the left and function of $t$ on the right, for them to be equal, they need to be equal to some constant (the constant is chosen negative to allow for a wave-like solutions):
$$\frac{X''}{X}+B(x)\frac{X'}{X}=-\mu^2$$
$$\frac{1}{c^2} \frac{T''}{T}=-\mu^2$$
The general solution for $T$ is simple:
$$T(t)=C_1 e^{i \mu c t}+C_2 e^{-i \mu c t}$$
Where $C_1$ and $C_2$ are some constants.
The second equation is an ODE and can be solved numerically for general $B(x)$:
$$X''+B(x)X'+\mu^2 X=0 \tag{4}$$