numerical solution of continuity equation, implicit scheme, staggered grid

Hi!

I'm trying to implement an implicit scheme for the continuity equation.

The scheme is the following:

http://img28.imageshack.us/img28/319...11130at003.png

With $\displaystyle \rho$ being the density, $\displaystyle \alpha$ is a weighing constant. d is a parameter that relates the grid spacing to the velocity of flow. j is the spatial grid and n is the time grid.

The problem with this scheme is the fact that it is implicit and effectively I have no idea how to successfully implement it. I tried by doing the following:

Assume have initial time n=0 spatial (j) grid full. Also assume I know spatial boundaries (j=-1/2 and j=max) at time= n+1.

Then I set j=1/2 equal to a variable x.

Next rearrange so that have j=3/2 in terms of variable x and boundary.

Do likewise for all the (half)grid points so that they can be written in terms of x. then when I reach the opposite boundary of grid (j=max) solve for x. So then I have j=1/2 and hence can substitute x into all the other equations to fill the grid points with data.

This can't work since:

When writing expression for any grid point in terms of x, i divide by $\displaystyle 0.5d(1-\alpha)$. Which is roughly 0.5. So going through all the grid points I end up with 0.5^400 (400 grid points) in the denominator which is bound to kill the calculation.

Can anyone tell me how to deal with this implicit scheme?