function [DX,DY,Lap] = diffOps_sol(h,N,bsq)
% DIFFOPS: construct the operators corresponding to computing partial
% derivatives and a Laplacian. The Laplacian is augmented with Dirichlet
% boundary conditions. 
%
% INPUTS: h: interval between two consecutive points
%         N: the grid is made of (N+1)^2 points
%         bsq: indices of the boundary points
%
% OUTPUTS: DX: discretization of \partial_x
%          DY; discretization of \partial_y
%          Lap: discretized Laplacian


% matrix of the 1-dimensional operator d/dx
vecD = [0 1 zeros(1,N-1)]/(2*h);
D = toeplitz(-vecD,vecD);

% correction of the boundary value
D(1, 1:3) = 1/2/h*[-3 4 -1]; D(N+1, N-1:N+1) = 1/2/h*[1 -4 3];

% sparsify
D = sparse(D); I = speye(N+1);

% operators d/dx and d/dy on a square grid
DX = kron(D,I);
DY = kron(I,D);

% Laplacian with Dirichlet boundary conditions
D2 = D*D; Lap = kron(I,D2) + kron(D2,I);
Lap(bsq,:) = zeros(4*N,(N+1)^2); Lap(bsq,bsq) = eye(4*N);

