% the problem is: 
%   find u:[0,1]x[0,1] -> R tha satisfies
%   -(D_xx + D_yy) u = f
%   u = 0 on the boundary of [0,1]x[0,1]
% The partial derivatives are discretized on a uniform (n-1)x(n-1) grid
% with gridsize h=1/n.
% b=x*(1-x)*h^2 
function [A, b]=get2Dbdvp(n)
h=1/n;
row=zeros(5*(n-1)^2,1);
col=zeros(5*(n-1)^2,1);
val=zeros(5*(n-1)^2,1);

% counter
ct=0;
for i=1:n-1
    for j=1:n-1
        % middle entry
        ct=ct+1;
        row(ct) = i+(j-1)*(n-1);
        col(ct) = i+(j-1)*(n-1);
        val(ct) = 4;
        % entry coresponding to the left value
        if i>1
            ct=ct+1;
            row(ct) = i+(j-1)*(n-1);
            col(ct) = i-1+(j-1)*(n-1);
            val(ct) = -1;
        end
        % entry coresponding to the right value
        if i<n-1
            ct=ct+1;
            row(ct) = i+(j-1)*(n-1);
            col(ct) = i+1+(j-1)*(n-1);
            val(ct) = -1;
        end
        % entry coresponding to the value below
        if j>1
            ct=ct+1;
            row(ct) = i+(j-1)*(n-1);
            col(ct) = i+(j-2)*(n-1);
            val(ct) = -1;
        end
        % entry coresponding to the left value
        if j<n-1
            ct=ct+1;
            row(ct) = i+(j-1)*(n-1);
            col(ct) = i+j*(n-1);
            val(ct) = -1;
        end
    end
end
row=row(1:ct);
col=col(1:ct);
val=val(1:ct);
b=zeros(1, (n-1)^2);
i=floor(n/2);
b(i+(i-1)*(n-1))=h^2;
ct
A=sparse(row, col, val);

b=b';