Your numbers for the crust of the star look pretty reasonable.
Gravitational forces will matter - gravity is strong enough at the surface that you will have to use general relativity to calculate with - it is likely about a ten percent correction to Newtonian gravity on the surface, and it is stronger than Newtonian gravity would be. And nucleon-nucleon forces and also quark-quark forces will definitely matter, too.
You should know that deep in the star the density will, of course, rise even more. The densities will be much higher than in the crust, certainly.
You should also know that the equation of state of matter at very high baryon density, like in a neutron star, is an outstanding and unsolved theoretical problem. Essentially none of the approximation methods that exist are at all reliable at very much above the density of ordinary nuclei, which is below infinite nuclear matter density.
So what really has to be done is to turn the problem around and ask what do the observations of neutron stars teach us about the equation of state - the equation of state relates the pressure to the density, or if you like the energy to the pressure and from that we can calculate the speed of sound by:
[math]v_s^2 = \frac{dp}{d\rho}\vert_S[/math],
where the derivative is taken at constant entropy. This would also reproduce your formula, the Newton-Laplace formula, for an elastic medium, but deep in a neutron star it won't be a Young's modulus it will be an adiabatic compressibility that comes in instead.
Unfortunately, we don't have any measurement of the radius of a neutron star at all, and we have very few good measurements of the masses, so we don't have very strong constraints on the very high density equation of state.
There is an absolute limit to the speed of sound if relativity is correct, due to causality. So we know for sure that:
[math]v_s \leq c[/math]
But it is not clear that any ordinary form of matter can actually saturate that bound or even come that close to it, which is what you are asking, I think.
It is not at all beyond the realm of possibility that the baryon densities in the core of a neutron star may approach 8-10 times the baryon density of symmetric nuclear matter.
In this case elementary considerations very strongly suggest that there would be quark matter in the core of the star. However lattice QCD is not easily capable of handling finite baryon density.
But asymptotically, at very high densities, one can argue that the equation of state of quark matter would be expected to approach that for a non-interacting ultra-relativistic (massless) gas, due to asymptotic freedom of QCD. For an ideal relativistic gas of massless particles the velocity of sound is given by:
[math]v_s^2 = \frac{1}{3}c^2[/math]
(This follows from the tracelessness of the energy-momentum tensor, which implies [math]p=\frac{1}{3}\rho c^2[/math].)
If the particles are massive instead then the velocity of sound will be less than this.
Real lattice calculations at high temperature, but zero baryon density, show a sound speed in quark matter that is a bit below this value, more like [math]v_s^2 = 0.3 c^2[/math]. This is because there is more than one degree of freedom - there are quarks and gluons as well. But the relevance for the neutron star is not entirely apparent since the lattice only works at zero baryon density. See figure 4 of the following paper for the temperature dependence of the sound speed in QCD from lattice calculations.
It's very easy to think of non-relativistic theories that violate the bound, but these will not be valid at high densities.
A theory with conformal invariance will, I believe, always satisfy this bound.
So I think that to violate it in a relativistic theory you would need to have strong interactions and to also violate conformal invariance. But I don't know of any tractable relativistic theory that satisfies those conditions - so I don't know of a theory in which you can really calculate that will give you a faster speed of sound than the above.
Now you can just write down an equation of state that will do it easily: you could just say that [math]p=\rho c^2[/math], but I don't think that this will be a realistic approximation to the equation of state in the core of a neutron star. It will be much softer than this because at very short distances the "hard core" repulsion between baryons will break down as the baryons begin to overlap and the quarks and gluons begin to see each other directly.
It doesn't mean that it can't happen in a neutron star that the speed of sound approaches the speed of light.
It just means I don't know of a way for it to get quite that high without the matter being something other than the kind of matter that I would expect to be there in the core of the neutron star!