% Numerical approximation to diffusion equation using backward difference
%
% Input:
%        ifunc : the initial value function
%            a : the length of the rod
%            b : the time interval
%            n : n+2 is the number of spatial points of the mesh
%            m : the number of time points of the mesh
%            alpha2 : diffusivity
% Ouput:
%            p : the numerical solution of diffusion equation at the mesh points
%          x,t : the x-t mesh.
%
function [u, x, t] = myfd1difbak(alpha2, ifunc, a, b, n, m)

	h = a/(n+1);   % spatial spacing
	dt = b/(m-1);   % time spacing
	[x, t] = meshgrid(0:h:a, 0:dt:b);   % Uniform mesh, including boundary points.

    mu = alpha2*dt/h^2;

	u = zeros(n+2, m);
	idx = 2:n+1;

	d1 = ones(n-1, 1);
	d1 = -mu*[d1; 0];
	d3 =  ones(n-1, 1);
	d3 = -mu*[0; d3];
	d2 = (1 + 2*mu)*ones(n, 1);
	A = spdiags([d1 d2 d3], [-1 0 1], n, n);

	% Compute initial values
	u(idx, 1) = ifunc(x(1, idx));

    for it = 2:m
	    u(idx, it) = A\u(idx, it-1);
    end

end
