For the efficient representation of discrete functions and for the solution of PDEs, the sparse grid technique has been developed in the last years. It is mainly based on the finite element approach using a specific tensor product of 1D hierarchical basis ansatz functions. The resulting linear system can be solved efficiently by multilevel methods. While usual discretization methods require basically O(h-d) grid points, the sparse grid approach needs only O(h-1 log(h-1)d-1) grid points. Here, h denotes the mesh size employed and d denotes the dimension of the problem. The accuracy of the approximation, however, deteriorates pointwise and with respect to the L2- and Lmax-norm only slightly from O(h2) to O(h2 log(h-1)d-1) provided that the function to be represented is sufficiently smooth. In the non-smooth case, adaptive refinement can be applied straightforwardly and helps to maintain the advantage of the sparse grid approach over an adaptive conventional h-version of the finite element method. However, the setup and solution of the linear system is quite complicated, especially in the case of PDEs with non-constant coefficient functions. A sparse grid approach using the finite difference philosophy is in this respect much more simple and gives in some cases even a better preformance, i.e. a O(h2) accuracy without the log-terms. However, the resulting systems are not more symmetric in general. Furthermore, the solver employs the BiCG iteration and a multilevel preconditioner using so called prewavelets. It converges independently of the mesh size of the problem. We report on our method and its application to 2D and 3D PDEs of 2nd order with general coefficient functions, discuss the adaptive refinement approach and describe how we can deal with general complicated domains. Furthermore, we present the prewavelet multilevel preconditioner in detail.