That seems to solve the problem. The big shift is now gone.
Did you also test yours on WGS84 or only on a sphere?
You also have to know that proj.4 did not have the accurate implementation of the projection algorithm until a colleague of mine added it in version 4.9.1. But that only produces an error of a few meters and not the kind of shift we saw.
From our internal tests these are correct values calculated by pyproj.
input coordinate pairs
+proj=longlat +datum=WGS84 +no_defs
output coordinate pairs
+proj=aeqd +ellps=WGS84 +lat_0=8.5 +lon_0=21.5 +x_0=5621452.01998 +y_0=5990638.42298