function [X, Y, U] = dirichlet(h, Region, fn, fnborder, varargin)

% Задача Дирихле для уравнения Пуассона
% $\Delta u = fn(x, y)$, $u |_{\partial \Gamma} = fnborder$

    if nargin < 1 || isempty(h)
        h = 0.02;
    end;
    if nargin < 2 || isempty(Region)
        Region = 'L';
    end;
    if nargin < 3 || isempty(fn)
        fn = @fun;
    end;
    if nargin < 4 || isempty(fnborder)
        fnborder = @funborder;
    end;

    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); % нумерация узлов области
    B = numgridborder(R);   % узлы, лежащие на границе
    Region = R > 0; % шаблон из нулей и единиц для все узлов области
    Border = B > 0; % шаблон из нулей и единиц для узлов на границе
    region = R(Region); % номера всех узлов в определенном порядке
    border = R(Border); % номера узлов, лежащих на границе
    nregion = length(region);
    nborder = length(border);

    A = -delsq(R)/h^2;

    % правая часть уравнения
    f = zeros(nregion, 1);
    f(region) = feval(fn, X(Region), Y(Region), varargin{:});    

    % граничные условия
    A(border, :) = sparse(1:nborder, border, 1, nborder, nregion);
    f(border) = feval(fnborder, X(Border), Y(Border), varargin{:});  

    u = A\f;
    U(Region) = u(region);

    X = X(2 : n - 1, 2 : n - 1);
    Y = Y(2 : n - 1, 2 : n - 1);
    U = U(2 : n - 1, 2 : n - 1);

    clf
    shg
    surfc(X, Y, U);
    shading interp;

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

    function f = funborder(x, y)
        % f = x - y;
        % f = sin(pi*x) - sin(pi*y);
        f = 0;