function [x,t,U] = solveWaveEq(f,T,u0,du0,g0,g1,h,tau,sigma)
% function [x,t,U] = solveWaveEq(f,T,u0,du0,g0,g1,h,tau,sigma)
% Solving the wave equation on rectangle [0,1]x[0,T] using Finite
% Differences Method - general scheme with weight sigma
% The equation:
%
%    d^2/dt^2 = d^2u/dx^2 + f(x,t)
%
% f      - annonymous function of variable x and t
% T      - upper bound for time
% u0,du0 - initial conditions (annonymous functions of variable x)
% g0,g1	 - boundary conditons (annonymous functions of variable t)
% h      - step in x
% tau    - step in time
% x,t    - knots in both directions including boundaries
% U      - approximate solution in the knots - matrix
%

N=round(1/h);
h=1/N;              % recounting of h
x=linspace(0,1,N+1);
xx=x(2:end-1);		% inner points
U0=u0(x);           % initial condition
UU0=U0(2:end-1);    % initial condition in inner points
UU0=UU0(:);
D0=du0(x);           % initial condition for the first derivative
DD0=D0(2:end-1);
DD0=DD0(:);
r=tau^2/h^2;
t=0;
tt=t;

dB=1+2^sigma*r;lB=-sigma*r;uB=-sigma*r;
dA0=-(1+2*sigma*r);lA0=sigma*r;uA0=sigma*r;
dA1=2*(1-(1-2*sigma)*r);lA1=(1-2*sigma)*r;uA1=(1-2*sigma)*r;

f0=tau^2*f(xx,0);
f0=f0(:);
B00=g0(0);          % boundary value in 0 for t=0
B0p=g0(tau);        % boundary value in 0 for t+tau
BD0=du0(0);         % boundary value in 0 for the derivative
B0m=B0p-2*tau*BD0;  % boundary value in 0 for t-tau counted from the derivative

B10=g1(0);          % boundary value in 1 for t=0
B1p=g1(tau);        % boundary value in 1 for t+tau
BD1=du0(1);         % boundary value in 1 for the derivative
B1m=B1p-2*tau*BD1;  % boundary value in 1 for t-tau counted from the derivative

f0(1)=f0(1)+sigma*r*B0m+sigma*r*B0p+(1-2*sigma)*r*B00;
f0(end)=f0(end)+sigma*r*B1m+sigma*r*B1p+(1-2*sigma)*r*B10;
R=mult3diag(dA1,lA1,uA1,UU0)-2*tau*mult3diag(dA0,lA0,uA0,DD0)+f0;
UU1=solve3diag(dB-dA0,lB-lA0,uB-uA0,R);

U=[B00 UU0' B10;B0p,UU1',B1p];
t=[t,tau];
tt=tau;

while tt<T
    tt=tt+tau; t=[t,tt];
    fk=tau*f(xx,tt); fk=fk(:);
    B0m=B00;B00=B0p;B0p=g0(tt); % boundary values for x=0
    B1m=B10;B10=B1p;B1p=g1(tt); % boundary values for x=1
    fk(1)=fk(1)+sigma*r*B0m+sigma*r*B0p+(1-2*sigma)*r*B00;
    fk(end)=fk(end)+sigma*r*B1m+sigma*r*B1p+(1-2*sigma)*r*B10;
    R=mult3diag(dA1,lA1,uA1,UU1)+mult3diag(dA0,lA0,uA0,UU0)+fk;
    UU2=solve3diag(dB,lB,uB,R);
    U=[U;B0p UU2' B1p];
    UU0=UU1;
    UU1=UU2;
end
    
end     % That's all, Folks!

