In principle, you cannot do this in a unique, accurate way without reconstructing something like the original DEM.
Let's see why, by doing the computation with a (difficult but realistic) example. It begins with a small portion of a 30 meter DEM, about 15 by 11 kilometers.
![Hillshaded DEM](https://i.stack.imgur.com/EHKmJ.png)
Contours have been computed at regular intervals: this is where you begin.
![Contours](https://i.stack.imgur.com/brWDv.png)
To perform the fill, we need to (1) identify the regions between contours and convert them to polygonal (or raster) features suitable for filling and (2) find the elevations of all contours bounding each region. Although this can be done entirely with vector operations, that would take some intermediate-level programming and is much trickier than it might seem. A raster solution, on the other hand, is relatively easy.
To obtain the regions between the contours, convert the contours to a grid format, construct the indicator grid of the NoData values, and extract its morphological components (aka "RegionGroup").
![Gaps](https://i.stack.imgur.com/X8959.png)
The purple regions are the gaps between contours.
![Regions](https://i.stack.imgur.com/lLieY.png)
Each region has been given a separate identifier and colored with a random color.
Returning to the contours, compute their Euclidean allocation grid. The value at each cell is the value of the contour closest to it.
![Euclidean allocation](https://i.stack.imgur.com/ZkOJc.png)
This image is zoomed into the top center of the previous ones to show detail. The contour lines are superimposed to show the relationship between the Euclidean allocation values and the nearest contours.
Finally, a zonal summary of the Euclidean allocation values, using the regions as zones, captures all the needed information: everything that can possibly be calculated without recreating a DEM is contained in this summary.
![Result](https://i.stack.imgur.com/kH20h.png)
The original hillshaded DEM shows in the margins beneath a choropleth map of the zonal minimum within each region (that is, the value of the lowest contour bounding each region). The contour lines are drawn on top for reference.
Success! The regions between the contour lines have all been filled according to the values of their bounding contours. So where are the alleged problems? Let's look at the zonal statistics table. This one has been sorted by decreasing area so we can focus first on the largest regions between the contours.
![Zonal summary table](https://i.stack.imgur.com/oFhmp.png)
See that range of zero in the fourth line? It is a large region which is bounded only by contours at 300 (meter) elevations. Let's highlight it in the map:
![Highlighted region](https://i.stack.imgur.com/dVBpd.png)
Most likely this region contains elevations between 200 and 300 meters, but the lower 200 meter contour does not appear anywhere on the map. (This particular region appears to be near the head of a valley which extends to the west off the map.) Because of this, the only contours with which it is in contact are at 300 meters. Now, how is the software supposed to be able to tell whether this is the upper contour or the lower one? The answer--which you have already determined visually--is that you must somehow extrapolate the elevation values downward from the peaks, across the 300 meter contour, and into lower levels. Doing that requires knowing how topography actually behaves and amounts to estimating elevation values in the between-contour regions. That means you have to recreate the DEM! (At least in a semi-quantitative manner: your extrapolation can be crude, provided it's in the right direction. For instance, if you happened also to have a layer of point elevations, with at least one point within each problem region, you could use that (via a spatial join) to classify every region accurately.) Without this additional calculation of elevations, it is not possible to determine whether this region should be classified as being between 200-300 meters or 300-400 meters.
How common is this problem? It occurs along all the edges as well as within any peak and any valley:
![Problem regions](https://i.stack.imgur.com/IGzoH.png)
In other words, it will always occur somewhere in a contour map.
I suppose you could go through and manually assign values to the problem regions according to your visual assessment of the map, but in most cases it would be easier, faster, and more reliable simply to create a DEM--even a crude one will do--and use that to classify and symbolize the regions.
All the calculations performed to make these illustrations can be carried out exactly as described, step-for-step, in QGIS/GRASS, ArcGIS/Spatial Analyst, and most other raster-based GISes. For details, see the help pages for your favorite GIS.
Best Answer
You don't want to "convert" vector contours to raster format: that would merely create a raster whose values designate which contour crosses each cell (and ought to have null values in cells not crossed by contours).
To convert a set of contours to a DEM you need to interpolate the values from the contours into the regions between and around the contours. According to the current GRASS manual, the proper procedure is r.surf.contour. It's a crude algorithm: using a vector contour dataset that has been converted to raster (as described in the question), it linearly interpolates between the nearest points found on the two contours immediately surrounding each cell. This process lops the tops off all peaks, fills lakes, flattens ridges and valleys, and tends to create partial step-like surfaces. (Think of an irregular staircase at the beach onto which a lot of sand has been blown: the edges of the steps are the contour lines and the sand fills the spaces between them.) Wherever the contour lines are closely and evenly spaced this will do a good job, but the results elsewhere may look less than satisfactory.
Other options
Better ways to interpolate directly between vector contour lines exist, such as the related technique I described at http://forums.esri.com/Thread.asp?c=3&f=39&t=102418&mc=20#288514. It uses multiple transects through each grid point to select neighboring contour locations and applies polynomial, rather than linear, interpolation to their values.
As Michael Miles-Stimson answered in a comment, you can apply techniques to triangulate the contour nodes and interpolate across the triangles. These sometimes work, but they will fail to capture peaks and valleys correctly (unless you have separate data to locate them, in which case triangulation is a singularly useful method) and their results will depend on the density of nodes used to represent the contours.
If you have access to ESRI software, use Topo To Raster. This automatically uses a combination of methods (including triangulation and some hydrological techniques) to produce interpolations that are physically realistic.
Another option is to use kriging with variable support--but that's an extremely specialized procedure that I believe is available only in expensive proprietary software used in the resource discovery and extraction industries (oil and mining).