Add fd stabilization

fd/fd_transport_stabilized.m

+function A = fd_transport_stabilized(epsilon, beta, T, n, timesteps, varargin)
+
+tau = T/timesteps;
+h = 1/(n-1);
+
+% Peclet number
+P = abs(beta*h/(2*epsilon));
+
+if (P>1)
+    B = - epsilon*(1+P)*fd_laplace(n) - beta*fd_gradient(n);
+else
+    B = - epsilon*fd_laplace(n) - beta*fd_gradient(n);
+end
+
+B = speye(n,n) + tau*B;
+
+if length(varargin)>0
+    bc_type=varargin{1};
+else
+    bc_type='dirichlet';
+end
+
+if strcmpi(bc_type,'dirichlet')
+    % Dirichlet boundary condition
+    B(1,:)=0;
+    B(1,1)=1;
+    B(end,:)=0;
+    B(end,end)=1;
+end
+
+if strcmpi(bc_type,'neumann')
+    % Neumann boundary condition
+    B(1,:)=0;
+    B(1,1)=-1/h;
+    B(1,2)=1/h;
+    B(end,:)=0;
+    B(end,end)=1/h;
+    B(end,end-1)=-1/h;
+end
+
+u0 = @(x) (x>=.4).*(x<=.5);
+
+xx = linspace(0,1,n)';
+u=xx*0;
+u(:)=u0(xx);
+
+u(1) = 0;
+u(end) = 0;
+
+
+plot(xx,u);
+a=axis();
+
+for k=1:timesteps
+    pause(.1);
+
+    rhs = u;
+    
+    if strcmpi(bc_type,'neumann')
+        % Neumann boundary condition
+        rhs(1) = 0;
+        rhs(end) = 0;
+    end
+
+    u = B\rhs;
+    plot(xx,u);
+    axis(a);
+end
+
+