function wave(h, tau, Region, fn, fnborder, fninitial, fndinitial, explicit, varargin)

% Начально-краевая задача для волнового уравнения разностным методом
% $\frac{\partial^2 u}{\partial t^2} = \Delta u + fn(t, x, y)$, $(x,y)\in \Gamma$
% $u(t, x, y) = fnborder(t, x, y)$, $(x,y)\in \partial\Gamma$
% $u(0, x, y) = fninitial(x, y)$
% $\frac{\partial u(0, x, y)}{\partial t} = fndinitial$
% $h$ --- шаг по пространственным координатам $x$, $y$
% $tau$ --- шаг по $t$
% explicit = 1 --- явная схема
% explicit = 0 --- неявная схема
% $Region$ --- форма области (см. $numgrid$)


    if nargin < 1 || isempty(h)
        h = 0.05;
    end;
    if nargin < 2 || isempty(tau)
        tau = 0.01;
    end;
    if nargin < 3 || isempty(Region)
        Region = 'L';
    end;
    if nargin < 4 || isempty(fn)
        fn = @fun;
    end;
    if nargin < 5 || isempty(fnborder)
        fnborder = @funborder;
    end;
    if nargin < 6 || isempty(fninitial)
        fninitial = @funinitial;
    end;
    if nargin < 7 || isempty(fndinitial)
        fndinitial = @fundinitial;
    end;
    if nargin < 8 || isempty(explicit)
        explicit = 1;
    end;

    stop = uicontrol('style', 'togglebutton', 'string', 'stop');

    x0 = 0;
    x1 = 1;
    n = (x1 - x0)/h + 3;
    [X, Y] = meshgrid(linspace(x0 - h, x1 + h, n));
    U = zeros(n); 
    U(:) = NaN;

    R = numgrid(Region, n);  % region grid 
    B = numgridborder(R);  % border grid
    A = -delsq(R);        % разностный оператор Лапласа

    Region = R > 0;      % шаблон из нулей и единиц
    region = R(Region);  % номера элементов, содержащих 1
    Border = B > 0;
    border = B(Border);

    len = max(max(R));

    u_cur = zeros(len, 1);
    dudt = zeros(len, 1);

    u_cur(region) = feval(fninitial, X(Region), Y(Region)); % U - матрица, u - вектор
    dudt(region) = feval(fndinitial, X(Region), Y(Region));
    u_prev = u_cur - tau*dudt;    % для начальных условий

    shg;
    graph = surf(X, Y, U);
    axis([x0 x1 x0 x1 -30 30]);

    t = 0;

    while ~get(stop, 'value')
      U(Region) = u_cur(region);
      set(graph, 'ZData', U);
      drawnow;

      t = t + tau;

      f = feval(fn, X(Region), Y(Region), t); 

      if explicit
         u_new = 2*u_cur - u_prev + tau^2/h^2 * A * u_cur + tau^2*f;
      else
         u_new = (speye(size(A)) - tau^2/h^2 * A ) \ (2*u_cur - u_prev + tau^2*f);
      end;
      u_prev = u_cur;
      u_cur = u_new;

      u_cur(border) = feval(fnborder, X(Border), Y(Border), t);
    end;

function f = fun(x, y, t)
    f = 500*(x - y); 

function f = funborder(x, y, t)
    f = 5*(x - y)*sin(2*pi*t);

function f = funinitial(x, y)
    f = 0;

function f = fundinitial(x, y)
    f = 0;