Objective To simulate the growth of axons in the developing nervous system numerically. Methods A set of mixed parabolic-gradient nonlinear equations was established based upon the principles of the developing neural biology, which was computed by the ADI scheme and the modified Euler' s method. Results (1) In the first case that one target was uniformly surrounded by ten axons in a square area, the pathway of the growing axons we simulated coincided with those obtained by the formers; (2) in the second case that a group of two-by-two-placed targets was surrounded by the same axons as in (1), the results correctly reflected the chemotaxis of the growing axons; (3) in the final case that a stream of targets lined near the upper edge of the square area and a stream of axons lined near the lower edge, the results approximately displayed the bundling and debonding courses of the growing axons. Conclusion The mathematical model and the numerical method presented in this paper can simulate the main phenomenon observed in the related experiments.