% generate the standard finite difference matrix with Dirichlet b.c.
A = sysmatrix_poisson;

% construct the Psi matrix (correction terms)
Psi_D = corrections(IX_D);
Psi_N = corrections(IX_N);

PSI_D = Psi_D.xx + Psi_D.yy;
PSI_N = Psi_N.xx + Psi_N.yy;

P = [PSI_D , PSI_N];

N = GRID.nx*GRID.ny;

% construct the D-matrix and the vector Ft (jump approximation
% with the one-sided least squares)
% boundary conditions go with '-' sign, as the third argument
% requires the jump.
[D_D,Ft_D] = d_matrix(IX_D,-uD,[],M,'Dirichlet'); 
[D_N,Ft_N] = d_matrix(IX_N,[],-uN,M,'Neumann');

D  = [D_D ; D_N ];
Ft = [Ft_D; Ft_N];

IN = length(IX.type);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the EJIIM system
EJIIM_MATR = [A P;D speye(6*IN,6*IN)];
EJIIM_F    = [F;Ft];

PSI = P;    

% solve the EJIIM system
SOL = EJIIM_MATR\EJIIM_F;

% separate solution and jump components
U = SOL(1:N);
J = SOL(N+1:end);