function [u, X, T, M, N] = ...
    WaveEquation1D(c, u0, v0, fL, fR, boundary, t, x, lambda, h)
% Równanie 
%     u_tt = c^2*u_xx 
% z warunkami pocz¹tkowymi 
%     u(t0,x) = u0(x)
%     u_t(t0,x) = v0(x)
% i warunkami brzegowymi fL(t) oraz fR(t) typu Dirichleta lub Neumanna na,
% odpowiednio, lewym i prawym brzegu. Parametry: 
% t = [t0, t1] - przedzia³ czasowy
% x = [x0, x1] - przedzia³ przestrzenny
% h - krok przestrzenny, lambda = h/k, gdzie k - krok czasowy.
% Zmienna boundary przyjmuje jedn¹ z dwóch wartoœci:
% (1) 'Dirichlet'
% (2) 'Neumann'
% w zale¿noœci od wymaganego typu warunku brzegowego.

k = lambda*h;

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

% zainicjowanie macierzy rozwi¹zania:
u = zeros(M,N);

% uzupe³nienie warunku pocz¹tkowego dla n = 1
for i=1:M, u(i,1) = feval(u0,X(i)); end
if strcmp(boundary,'Dirichlet')     % punkty "widmowe" poza siatk¹
    uL = 2*feval(fL,T(1)) - u(1,1);
    uR = 2*feval(fR,T(1)) - u(M,1);
elseif strcmp(boundary,'Neumann')
    uL = u(1,1) - h*feval(fL,T(1));
    uR = u(M,1) + h*feval(fR,T(1));
end

% uzupe³nienie warunku pocz¹tkowego dla n = 2
for m=2:(M-1)
    u(m,2) = u(m,1) + k*feval(v0,X(m)) + ...
        c^2*lambda^2*(u(m+1,1)-2*u(m,1)+u(m-1,1))/2;
end
    % na brzegu korzystamy z punktów "widmowych" z poprzedniego kroku:
    u(1,2) = u(1,1) + k*feval(v0,X(1)) + c^2*lambda^2*(u(1+1,1)-2*u(1,1)+uL)/2;
    u(M,2) = u(M,1) + k*feval(v0,X(M)) + c^2*lambda^2*(uR-2*u(M,1)+u(M-1,1))/2;
if strcmp(boundary,'Dirichlet')     % punkty "widmowe" w aktualnym kroku:
    uL = 2*feval(fL,T(2)) - u(1,2);
    uR = 2*feval(fR,T(2)) - u(M,2);
elseif strcmp(boundary,'Neumann')
    uL = u(1,2) - h*feval(fL,T(2));
    uR = u(M,2) + h*feval(fR,T(2));
end

% g³ówna pêtla
for n = 3:N
    for m = 2:(M-1)
        u(m,n) = 2*u(m,n-1) - u(m,n-2) + ...
            c^2*lambda^2*(u(m+1,n-1)-2*u(m,n-1)+u(m-1,n-1));
    end
    % na brzegu korzystamy z punktów "widmowych" z poprzedniego kroku:
    u(1,n) = 2*u(1,n-1) - u(1,n-2) + ...
        c^2*lambda^2*(u(1+1,n-1)-2*u(1,n-1)+uL);
    u(M,n) = 2*u(M,n-1) - u(M,n-2) + ...
        c^2*lambda^2*(uR-2*u(M,n-1)+u(M-1,n-1));
    
    if strcmp(boundary,'Dirichlet') % punkty "widmowe" w aktualnym kroku:
        uL = 2*feval(fL,T(n)) - u(1,n);
        uR = 2*feval(fR,T(n)) - u(M,n);
    elseif strcmp(boundary,'Neumann')
        uL = u(1,n) - h*feval(fL,T(n));
        uR = u(M,n) + h*feval(fR,T(n));
    end
    
end