Hi everyone: I am rookie in classical physics and first-time PF user so please forgive me if I am making mistakes here. My current project needs some guidance from physics and I am describing the problem, my understanding and question as below. I have an independent electrostatic system with positive and negative charges distributed on a 2D plane. The plane is discretized into small sub-rectangles to calculate the charge density and get the field and potential distribution using Poisson equation, so is the total potential energy N. Location of each charge c (suppose to be x(c)) determines N. Now I need to calculate the Hessian matrix (2nd-order partial differentiation) of N, H(N), with respect to each x(c). To me, as ∂^2 N /∂ x(c)^2 = -q(c)D(c), we get the diagonal elements of H(N), where q(c) is the quantity of charge and D(c) is the local density; the off-diagonal elements need to be calculated using Coulomb's law. Please correct me if there is any mistake in the above, e.g., miscalculation of H(N), as well as others. If all is correct, my question is (1) what kind of properties would H(N) have (I assume it would be symmetric, but would it be positive definite also)? (2) if we scale up the system with more charges, is there a fast way to approximate H(N)? Thanks a lot in advance.