This work presents an approach based on Neural Networks (NN) for constructing models of high speed bipolar transistors, to be used inside circuit simulators. Since recently physical modeling has been the prevailing approach to semiconductor modeling for circuit simulation, being Gummel-Poon the choice for bipolar junction transistors (BJT) at high frequencies. The knowledge of the underlying physical principles was sufficient for deriving effective models, where each parameter has a clear physical meaning. This situation is gradually changing, with the current size of bipolar junctions made available by the microelectronics manufacturing capabilities. As long as the ratio of surface to volume increases for the smallest devices, many boundary effects should be taken into account, leading to complex multidimensional analysis of the phenomena. Reliable models requires now a high development cost. Furthermore, precise physical models should now rely on a number of largely empirical fit parameters. The method here investigated is based on NN, which can be seen as a fitting techniques, with the advantage of output functions with unlimited degree of continuity. NN approaches has gained widespread attention over the past two decades, with a variety of applications, including dynamic systems identification, control and simulation. However, the traditional schemes used for dynamic NN, like tapped-delays feed-forward networks, feedback and self-recurrent networks, are not immediately usable for models to be embedded in circuit simulators. Simulators require a representation in terms of continuous-time equations, since time discretisation will be an internal process of the simulator itself. In the NN here adopted, nodes can be instances of differential equations. The whole network has two input nodes, the BJT voltages, and two output nodes, the BJT currents, and several internal layers populated with neurons. For the determination of the parameters a combination of training sets can be used, mixing measured responses of the device in DC conditions, small-signals in frequency domain, and large signals transient responses. A combination of global optimization followed by a quasi-Newton conjugate-gradient method has been used, for computing the internal parameters of the NN. Once the NN has been trained, the final set of parameters is used for generating a piece of code in Anacad's HDLA language, implementing the NN in recall mode. This will be the block modeling the BJT inside the ELDO circuit simulator.