function u = computeFDsolution(epsilon,b,c,f,bc,mesh,sigma_fct)
% computes the solution of FD with artificial diffusion to 
%
%      - epsilon*u'' + b*u' + c*u = f
% with boundary conditions 
%       bc(1,1)*u(0)+bc(1,2)*u'(0) = bc(1,3)
%       bc(2,1)*u(0)+bc(2,2)*u'(0) = bc(2,3)
%
% on a mesh given by mesh
% with artifical diffusion given by the function handle sigma

% boundary conditions 
alpha0 = bc(1,1);
beta0 = bc(1,2);
gamma0 = bc(1,3);
alpha1 = bc(2,1);
beta1 = bc(2,2);
gamma1 = bc(2,3);


% compute stiffness matrix
N=size(mesh,1)-1;
h = mesh(2:end)-mesh(1:(end-1));
hj = h(2:end);
hjm1 = h(1:(end-1));
hjav = 0.5*(hj+hjm1);

Ajj = -1./(hjm1.*hjav);
Ajjp1 = 2./(hj.*hjm1);
Ajjp2 = -1./(hj.*hjav);
A = spdiags([Ajj,Ajjp1,Ajjp2],[0,1,2],N-1,N+1);

Bjjp1 = -b(mesh(2:(end-1)))./hj;
Bjjp2 = -Bjjp1;
Bplus = spdiags([Bjjp1,Bjjp2],[1,2],N-1,N+1);

Bjj = -b(mesh(2:(end-1)))./hjm1;
Bjjp1 = -Bjj;
Bminus = spdiags([Bjj,Bjjp1],[0,1],N-1,N+1);

Bjj = -(hj./hjm1).*b(mesh(2:(end-1)))./(hj+hjm1);
Bjjp1 = -((hjm1./hj)-(hj./hjm1)).*b(mesh(2:(end-1)))./(hj+hjm1);
Bjjp2 = (hjm1./hj).*b(mesh(2:(end-1)))./(hjm1+hj);
B0 = spdiags([Bjj,Bjjp1,Bjjp2],[0,1,2],N-1,N+1);

C = spdiags(c(mesh(2:(end-1))),[1],N-1,N+1);

% change for different method
B = Bminus;

% add artificial diffusion
q = b(mesh(2:(end-1))).*hj/(2*epsilon);

sigma = sigma_fct(q);

S = epsilon*spdiags(sigma,0,N-1,N-1)*A+B+C;

F = f(mesh(2:(end-1)));


% include boundary conditions
% alpha0 u(0) + beta0 u'(0) = gamma0
% alpha1 u(1) + beta1 u'(1) = gamma1

Sbc = [alpha0 - beta0/h(1), beta0/h(1), sparse(1,N-1);
       S;
       sparse(1,N-1), -beta1/h(end), alpha1 + beta1/h(end)];

Fbc = [gamma0;F;gamma1];

% solve linear system
u = Sbc\Fbc;

end


