function [u, X, T, M, N] = DiffusionEquation(D, u0, fL, fR, t, x, mu, h)
% Równanie 
%     u_t = (D(x)*u_x)_x
% z warunkiem pocz¹tkowym
%     u(t0,x) = u0(x)
% i warunkami brzegowymi typu Dirichleta
%     u(t,x0) = fL(t)
%     u(t,x1) = fR(t).
% Parametry:
% t = [t0, t1] - przedzia³ czasowy
% x = [x0, x1] - przedzia³ przestrzenny
% h - krok przestrzeny, mu = k/(h^2), gdzie k - krok czasowy.
% NAG£ÓWEK DO ZMIANY (WAR. BRZEGOWE)

k = mu*(h^2);

% punkty siatki (siatka przestrzenna - centralna)
X = (x(1)+h/2):h:(x(2)-h/2);    M = length(X);
T = t(1):k:t(2);                N = length(T);

% zainicjowanie macierzy rozwi¹zania (uwzglêdniamy poza siatk¹)
u = zeros(M+2,N);

% uzupe³nienie warunku pocz¹tkowego dla n = 1
% LINIE 27-28 WYMAGAJ¥ ZMIAN
u(2:(M+1),1) = u0(X);
u(1,1) = 2*feval(fL,T(1)) - u(2,1);
u(M+2,1) = 2*feval(fR,T(1)) - u(M+1,1);

% wspó³czynniki uk³adu równañ dla punktów wewnêtrznych
alpha = -mu/2*feval(D,X-h/2);
beta = 1 + mu*(feval(D,X+h/2)+feval(D,X-h/2))/2;
gamma = -mu/2*feval(D,X+h/2);

% LINIE 36-38 WYMAGAJ¥ MODYFIKACJI (RÓ¯NE WAR. BRZEGOWE)
a = [0, alpha, 0.5]; % dolna przek¹tna macierzy uk³adu
b = [0.5, beta, 0.5]; % g³ówna przek¹ta macierzy uk³adu
c = [0.5, gamma, 0]; % górna przek¹tna macierzy uk³adu

delta = zeros(1,M); 

for n = 1:N-1
    % wspó³czynniki wektora "prawej strony" uk³adu równañ
    for m = 1:M
        delta(m) = u(m,n)*(mu/2*feval(D,X(m)-h/2)) + ...
            u(m+1,n)*(1-mu/2*feval(D,X(m)+h/2)-mu/2*feval(D,X(m)-h/2)) + ...
            u(m+2,n)*(mu/2*feval(D,X(m)+h/2));
    end
    % LINIA 50 WYMAGA ZMIAN
    d = [feval(fL,T(n+1)), delta, feval(fR,T(n+1))];
    u(:,n+1) = ThomasSolve(a,b,c,d);
end