Expression for electric potential around an electric dipole in Cartesian coordinates I'm trying to find the equation of the electric potential of a electric dipole in Cartesian coordinates. 
The equation should take the form of something looks like this:
$$z = x\exp(-x^2-y^2)$$

Where x and y are Cartesian coordinates and $z$ is the electric potential, and $-\nabla z$ is the electric field (quivers in the plot). 
I know that equation is not true for dipole electric potential. It's just there to show how the plot of the real equation should look like.
 A: I take it that the object you want to model is a finite dipole, i.e. one with a nonzero separation $d$ between two finite charges $\pm q$.
The potential for this is trivial to construct out of its basic constituents: the configuration has two point charges a distance $d$ apart, which means that the electric potential will just be the sum of the potentials produced by each of the charges. Setting the inter-charge axis along the $z$ axis, the potential reads
$$
V(x,y,z) = \frac{q}{\sqrt{x^2+y^2 +(z-d/2)^2}} - \frac{q}{\sqrt{x^2+y^2 +(z+d/2)^2}}.
$$

Here it's worth remarking on just how different this is to the expression you wrote down, particularly as regards its decay at large distances from the origin, which is algebraic (specifically, $V(\mathbf r) \sim qd/r^2$ for $r\gg d$) instead of exponential. Generally speaking, electrostatic potentials which decay exponentially as functions of $r$ are so uncommon that if you find such a potential and you don't know exactly what you're doing, you're basically guaranteed to be doing something wrong.
Now, I can't go as far as saying that the potential you wrote down is unphysical, because for every smooth electrostatic potential $V(\mathbf r)$ there is always a charge density $\rho(\mathbf r) := -\varepsilon_0 \nabla^2 V(\mathbf r)$ that induces it. However, it is important to call out some extremely unlikely features of the trial potential you've written down, which relate to the charge density which generates it (via the Laplacian). That charge density looks like this:

At first glance this looks reasonable, with two big lumps of opposite charge and some 'ringing' around them, but that appearance is deceptive:


*

*That 'ringing', i.e. the shallow regions of charge opposite to the central lumps at $r>1$, are much wider than the central lumps, and because of that, they have a significant fraction of the charge of central lumps themselves, but with opposite sign and at a larger distance from the origin (so they contribute more to the dipole moment).

*Moreover, because of that, the dipole moment of the charge distribution is zero.
(And, even worse, all the multipole moments of the charge distribution vanish. This is easy to see directly from the exponential decay of the potential: if any of the multipole moments were nonzero, it would induce an algebraic decay in the potential that's much slower than the exponential.)
I imagine none of these features are things that you initially contemplated. Unless you were actively seeking out a potential with all of these properties, an exponential decay in the potential is simply not appropriate.
A: If you want just a "D plot of the potential in the x-y plane which contains the changes then a simple way of getting a 3D representation and a contour plot is to use WolframAlpha.  
If you have equal magnitude positive and negative charges at $(-1,0)$ and $(+1,0)$ (with potentials inversally proportional to the distance from the charge) then the resulting potential (the algebraic sum of the potentials due to the individual charges) is of the form $$V(x,y) = \dfrac {-1}{\sqrt {(x-1)^2 + y^2}} + \dfrac {1}{\sqrt {(x+1)^2 + y^2}}$$
WolframAlpha input is z= -1/sqrt ((x-1)^2 + y^2) + 1/sqrt ((x+1)^2 + y^2) 
WolframAlpha you will produce the 3D representation with the potential in the z-direction and also a contour plot.  

Other applications will give you more control over the output.
