I am a newbie in finite difference methods, so I apologize in advance if the question is trivial.

I am trying to solve the advenction equation, i.e.
$\frac{\partial \phi(x,t)}{\partial t} + v \frac{\partial \phi(x,t)}{\partial x}=0$, where $v$ is a velocity, with the finite difference method. The system is of course also solvable analytically, and it turns out that every function of the kind $f(x,t) = f(x-vt)$ is a solution. Therefore, given an arbitrary initial condition $\phi(x,0)$, the "shape" of $\phi(x,0)$ will move towards positive x (for $v>0$) or negative x (for $v<0$).

I first formulated the problem with the backward formula for the spatial derivative, 

$\phi(x\_i,t\_{j+1}) = \phi(x\_i,t\_j) - v\cdot\Delta T\cdot(\phi(x\_i,t\_j)-\phi(x\_{i-1},t\_j))/\Delta x$.

The initial condition is given by a Gaussian packet centered at the origin,

$\phi(x\_i,t\_1) = exp(-(x\_i/2\sigma)^2)$.

With $v>0$ the gaussian packet does indeed move towards right with increasing times, and maintaining always the same shape. However, if I set $v<0$, I obtain large instabilities, with an exponentially increasing amplitude. If I reformulate the problem in terms of the forward formula for the spatial derivative 

$\phi(x\_i,t\_{j+1}) = \phi(x\_i,t\_j) - v\cdot\Delta T\cdot(\phi(x\_{i+1},t\_j)-\phi(x\_{i},t\_j))/\Delta x$

I obtain the opposite behavior: with a negative $v$, the packet propagates nicely towards left, but a positive $v$ leads to a large divergence of the system.

Finally, I tried to use the centered formula for the spatial derivative, 
$\phi(x\_i,t\_{j+1}) = \phi(x\_i,t\_j) - v\cdot\Delta T\cdot(\phi(x\_{i+1},t\_j)-\phi(x\_{i-1},t\_j))/(2\Delta x)$

and in this case I obtain divergences for both signs of $v$.

Is there a particular reason for this kind of behaviour? What should be the "most correct" way to implement this problem?

I attach here a short Matlab script that I made, and that reproduces what I described.

```
            L = 1; %Length
            v = 1; %Speed
            i=sqrt(-1);
            DeltaT = ((L/abs(v))/1000)
            DeltaX = 1*abs(v)*DeltaT
            x=-L:DeltaX:L;
            t = 0:DeltaT:(L/abs(v));
            sigma = L/50;
            Tbreak = 1e-5; %Parameter for animation

            %% v = 1, backward formula
            v=1;
            phi = zeros(length(x),length(t));phi(:,1) = exp(-((x)./(2*sigma)).^2);
            for n=1:(length(t)-1)
                for k=2:(length(x)-1)
                    phi(k,n+1)=phi(k,n) - v*DeltaT*(phi(k,n)-phi(k-1,n))/(DeltaX);
                end
            end
            figure(1);
            for n=1:(length(t))
            plot(x,phi(:,n))
            pause(Tbreak)
            end

            %% v = -1, backward formula
            v=-1;
            phi = zeros(length(x),length(t));phi(:,1) = exp(-((x)./(2*sigma)).^2);
            for n=1:(length(t)-1)
                for k=2:(length(x)-1)
                    phi(k,n+1)=phi(k,n) - v*DeltaT*(phi(k,n)-phi(k-1,n))/(DeltaX);
                end
            end
            figure(1);
            for n=1:(length(t))
            plot(x,phi(:,n))
            pause(Tbreak)
            end

            %% v = 1, forward formula
            v=1;
            phi = zeros(length(x),length(t));phi(:,1) = exp(-((x)./(2*sigma)).^2);
            for n=1:(length(t)-1)
                for k=2:(length(x)-1)
                    phi(k,n+1)=phi(k,n) - v*DeltaT*(phi(k+1,n)-phi(k,n))/(DeltaX);
                end
            end
            figure(1);
            for n=1:(length(t))
            plot(x,phi(:,n))
            pause(Tbreak)
            end

            %% v = -1, forward formula
            v=-1'
            phi = zeros(length(x),length(t));phi(:,1) = exp(-((x)./(2*sigma)).^2);
            for n=1:(length(t)-1)
                for k=2:(length(x)-1)
                    phi(k,n+1)=phi(k,n) - v*DeltaT*(phi(k+1,n)-phi(k,n))/(DeltaX);
                end
            end
            figure(1);
            for n=1:(length(t))
            plot(x,phi(:,n))
            pause(Tbreak)
            end

```