When solving the Poisson equation by means of centered finite difference, one usually ends up with a formulation of the type
$$ Au = f, $$
where $A$ encodes the finite difference coefficients and the boundary conditions.
Now, assume I want to solve the Poisson equation on the unitary square. I know that there exists a formula linking the 1D matrices and the overall 2D matrix, and such formula reads as
\begin{equation} A^{(2)} = I_{n_y} \otimes A_{x} + A_{y} \otimes I_{n_x}, \quad(1) \end{equation}
where $A_x$ (resp. $A_y$) is the matrix relative to the 1D discretization over the x-axis (resp. y-axis), $I$ is the identity matrix, and $n_x$ (resp. $n_y$) is the number of points over the x-axis (resp. y-axis).
I am interested in a general formulation where one can set Dirichlet or Neumann boundary conditions on any of the four sides of the squared domain, but I am not sure if it suffices to set properly $A_x$, $A_y$ and then rely on (1).
For example, assume I want to set Neumann on the left and on the top side, and Dirichlet on the remaining ones. I wonder if it sufficies to set a 1D problem (x-axis) with Neumann condition on the left point, Dirichlet on the right one, and another 1D problem (y-axis) with Neumann on the right and Dirichlet on the left, and then use (1) to reconstruct the 2D problem matrix.
I tried to figure it out mathematically through (1) using the Einstein notation, but I cannot understand, after a certain point, how to evaluate the boundary coefficient. As well, I tried to use third party codes (Python) to compare experimentally the obtained matrices, but I can't find a solver for 2D squares with general boundary conditions (I prefer to use "untouched" codes already validated to avoid error replication).