Hi! I'm trying to implement an implicit scheme for the continuity equation. The scheme is the following: http://img28.imageshack.us/img28/3196/screenshot20111130at003.png [Broken] With [TEX]\rho[/TEX] being the density, [TEX]\alpha[/TEX] 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 [tex]0.5d(1-\alpha)[/tex]. 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?