function [x,u] = diff1d_implicit(a,b,u0,K,mua,dt,nt)
% solve the 1D diffusion equation
% - d/dx (K dy/dx) + mua y(x) = dy/dt
% K,mua are constants
% u0 is initial function
% dt is time step
% nt is number of iterations
%
s = size(u0);
n = s(2)
h = (b-a)/(n-1);   % grid spacing
x = [a : h : b];  % n samples
A = x; B = x; C = x; S = x;
u = zeros(n,nt);
u(:,1) = u0';
for i = 1:n
% create the "full" diagonals
A(i) = -K*dt/(h*h);
B(i) = 1 + 2*dt*K/(h*h) + dt*mua;
C(i) = -K*dt/(h*h);

end

% make a sparse matrix
B = [A',B',C'];
d = [-1,0,1];
%S
Dmat = spdiags(B,d,n,n);
%full(K);
for k = 2:nt
    u(:,k) = Dmat\u(:,k-1);
end