function U=Solve2dR(f,ai,bi,gl,gr,p1,p2)
% SOLVE2DR solves 2d Poisson problem using Robin conditions
%   u=Solve2dR(f,ai,bi,gl,gr,p1,p2) solves the two dimensional
%   Poisson equation Delta u=f on the domain Omega=(ai*h,bi*h)x(0,1) with
%   Robin boundary conditions (dn+p1)u=gl at x=ai*h and (dn+p2)u=gr at
%   x=bi*h and u=0 at y=0 and y=1 using a finite difference
%   approximation with interior grid points (bi-ai) times length(gl)
%   using the same mesh size h=1/(length(gl)+1) in both x and y.

nx=bi-ai+1;                                  % number of point in x direction
ny=length(gl);                               % number of point in y direction
h=1/(ny+1);                                  % mesh size
ex=ones(nx,1);                               % construct 1d finite differences
Dxx=spdiags([-ex/h^2 2/h^2*ex -ex/h^2],[-1 0 1],nx,nx);
ey=ones(ny,1);
Dyy=spdiags([-ey/h^2 2/h^2*ey -ey/h^2],[-1 0 1],ny,ny);
A=kron(speye(size(Dxx)),Dyy)+kron(Dxx,speye(size(Dyy)));
A(1:ny,1:ny)=A(1:ny,1:ny)/2+p1/h*speye(ny);  % put Robin conditions
A(end-ny+1:end,end-ny+1:end)=A(end-ny+1:end,end-ny+1:end)/2+p2/h*speye(ny);
f(1:ny,1)=f(1:ny,1)/2+gl/h;                  % add boundary conditions into rhs
f(1:ny,end)=f(1:ny,end)/2+gr/h;
u=A\f(:);                                    % solve by sparse Gaussian el.
U=reshape(u,ny,nx);                          % put solution into matrix
