We describe an F-statistic search for continuous gravitational waves from galactic white-dwarf binaries in simulated LISA Data. Our search method employs a hierarchical template-grid based exploration of the parameter space. In the first stage, candidate sources are identified in searches using different simulated laser signal combinations (known as TDI variables). Since each source generates a primary maximum near its true Doppler parameters (intrinsic frequency and sky position) as well as numerous secondary maxima of the F-statistic in Doppler parameter space, a search for multiple sources needs to distinguish between true signals and secondary maxima associated with other, louder signals. Our method does this by applying a coincidence test to reject candidates which are not found at nearby parameter space positions in searches using each of the three TDI variables. For signals surviving the coincidence test, we perform a fully coherent search over a refined parameter grid to provide an accurate parameter estimation for the final candidates. Suitably tuned, the pipeline is able to extract 1989 true signals with only 5 false alarms. The use of the rigid adiabatic approximation allows recovery of signal parameters with errors comparable to statistical expectations, although there is still some systematic excess with respect to statistical errors expected from Gaussian noise. An experimental iterative pipeline with seven rounds of signal subtraction and re-analysis of the residuals allows us to increase the number of signals recovered to a total of 3419 with 29 false alarms.