function [x,y,U] = solveElipEq(p,q,f,g,R,mn)
% function [x,y,U] = solveElipEq(p,q,f,g,R,mn)
% Solving linear eliptical equation on rectangle using Finite Differences Method
% The equation:
%
%    -d(p(x,y)du/dx)/dx - d(p(x,y)du/dy)/dy + q(x,y)u = f(x,y)
%
% !!! Be careful of the sign '-' on the previous line. !!!
%
% p,q,f     	- annonymous functions of two variable
% g		- function for calculation the boundary values
% R=[a,b,c,d]   - the boudaries of the rectagle [a,b]x[c.d]
% mn=[m,n]	- number of division in both directions,
%		  one value may be used for m==n
%
% x,y  - knots in both directions including boundaries
% U    - approximate solution in the knots - matrix
%
% The points between knots are used for the approximation
% to obtain the symmetric matrix.
%

if length(mn)==1
    m=mn;n=mn;
else
    m=mn(1);n=mn(2);
end

% generating the knots
a=R(1);b=R(2);c=R(3);d=R(4);
Ix=[a,b];
Iy=[c,d];
x=linspace(a,b,m+1);
y=linspace(c,d,n+1);
xx=x(2:end-1);		% inner points
yy=y(2:end-1);
[XX,YY]=meshgrid(xx,yy);
hx=(b-a)/m;		% steps in both directions
hy=(d-c)/n;
x1=(x(1:end-1)+x(2:end))/2; % midpoints of the subintervals
y1=(y(1:end-1)+y(2:end))/2;
Q=q(XX,YY);
F=f(XX,YY);

U=zeros(n+1,m+1);   % matrix of the solution

% boundary conditions:
[Xp,Yp]=meshgrid(a,y);
U(:,1)=g(Xp,Yp);
[Xp,Yp]=meshgrid(b,y);
U(:,end)=g(Xp,Yp);
[Xp,Yp]=meshgrid(x,c);
U(1,:)=g(Xp,Yp);
[Xp,Yp]=meshgrid(x,d);
U(end,:)=g(Xp,Yp);

% right side of the system
F=f(XX,YY);
b=reshape(F',(m-1)*(n-1),1);

A=zeros((m-1)*(n-1),(m-1)*(n+1));   % matrix of the system of equations for inner points
                                    % The matrix is extended because of easier creation.
xL=x1(1:end-1);yL=yy(1);
[XL,YL]=meshgrid(xL,yL);
xLL=xx;yLL=y1(1);
[XLL,YLL]=meshgrid(xLL,yLL);
PL=diag(p(XL,YL))/(hx^2);                  % subdiagonal of the diagonal block
PLL=diag(p(XLL,YLL))/(hy^2);               % subdiagonal block
b(1:m-1)=b(1:m-1)+PLL*[U(1,2:end-1)]';     % boundary conditions for y=c are moved to the right side
for ii=1:n-1
    xR=x1(2:end);yR=yy(ii);
    [XR,YR]=meshgrid(xR,yR);
    xRR=xx;yRR=y1(ii+1);
    [XRR,YRR]=meshgrid(xRR,yRR);
    PR=diag(p(XR,YR))/(hx^2);              % superdiagonal of the diagonal block
    PRR=diag(p(XRR,YRR))/(hy^2);           % superdiagonal block
    ini=(ii-1)*(m-1);                      % initial index
    A(ini+1:ini+m-1,ini+1:ini+m-1)=-PLL;
    M=PL+PLL+PR+PRR;                       % diagonal of the diagonal block
    dL=diag(PL);                           % subdiagonal of the diagonal block
    dR=diag(PR);                           % superdiagonal of the diagonal block
    Qd=diag(Q(ii,:));
    M=M-diag(dL(2:end),-1)-diag(dR(1:end-1),1)+Qd;
    A(ini+1:ini+m-1,ini+(m-1)+1:ini+2*(m-1))=M;
    A(ini+1:ini+m-1,ini+2*(m-1)+1:ini+3*(m-1))=-PRR;
    b(ini+1)=b(ini+1)+dL(1)*U(ii+1,1);             % boundary conditions for left and right side of the rectangle
    b(ini+m-1)=b(ini+m-1)+dR(end)*U(ii+1,end); % are moved to the right side of the system
    PLL=PRR; % blocks are moved to the other side in the next step
    PL=PR;
end
b(ini+1:end)=b(ini+1:end)+PRR*[U(end,2:end-1)]'; % boundary conditions for y=d are moved to the right side

A=A(:,m:end-m+1);
U1=A\b;
UU=reshape(U1,m-1,n-1);
U(2:end-1,2:end-1)=UU';

end     % That's all, Folks!

