How to test the physical accuracy of a finite difference solution to Laplace equation? I just wrote a program that computes the potential field inside a Penning-Malmberg trap; my end goal is to simulate plasma inside this trap.
First I need to make sure my potential field across the grid is physically accurate. How can I do this? Which tests would you recommend me to perform on my grid?
I solved the field using a finite difference method; how can I determine what the appropriate grid size would be? Is a linear interpolation between grid points good enough to make good physical sense?
 A: 
How can I do this, which tests would you recommend me to perform on my grid?

When developing new code, the standard thing to do to ensure your code works as expected is to compare it to known solutions. For some mathematical models, there are closed form expressions (i.e., analytic solutions) and you can use these as a basis for testing. 
For instance, with the heat equation, you can initialize with a delta function and expect a Gaussian at later times (cf Wikipedia again). If we call the analytic solution $y(x,t)$ and the numerical solution $u(x,t,N)$ where $N$ represents the number of grid points, then the error (how far off your numerical model is from the analytic solution) is then
$$
E(t,N)=\left(\frac{1}{N}\sum_i\left[y(x_i,t)-u(x_i,t,N)\right]^2\right)^{1/2}\tag{1}
$$
(this is the $\ell^2$ norm, there are other choices). However, if you do not have an analytic solution, then you can compare your model to someone else's to the same effect. 

How can I determine what the appropriate grid size would be?

Or, for testing convergences and resolution, you can compare the difference of $N$ grid points versus $N/2$ and $2N$ grid points, as I mention in my answer to this question, using the same error condition in Equation (1). Note that you'll also run into issues of run-times with larger $N$, so the trade-off between convergence and solutions, so only testing can provide an answer to this question.
The ratio of the differences of the errors should be linear in $N$, so $E(N)/E(N/2)\sim N/(N/2)=2$. It could be different, however, for higher-order schemes.

Is a linear interpolation between grid points good enough to make good physical sense?

Finite difference methods assume constant values between grid points (i.e., the value of the potential/density/field/etc is grid-centered & approximately the same throughout). This of course is an approximation and can be mitigated by choosing a larger grid size. However, one can also assume some sort of interpolation (linear and quadratic are often chosen, but there's also "WENO" methods that are popular) of grid values, though this method is often referred to as finite volume, rather than finite difference.
Similar to the above, only testing can give you an answer as to whether something is "good" enough.
