Just to make sure we're all on the same page, I'll take a step back before taking two steps forward. The normalized wavefunctions for the $n = 2$ energy level are usually written as
\begin{gather}
\psi_{2,1,0}(r, \theta, \phi) = \left(\frac{1}{32\pi a_0^3}\right)^{1/2} \frac{r}{a_0} \mathrm{e}^{-r/2a_0} \cos(\theta), \\
\psi_{2,1,\pm1}(r, \theta, \phi) = \mp \left(\frac{1}{64\pi a_0^3}\right)^{1/2} \frac{r}{a_0} \mathrm{e}^{-r/2a_0} \sin(\theta) \mathrm{e}^{\pm\mathrm{i}\phi},
\end{gather}
where $a_0 = \hbar^2/m_\mathrm{e}e$ is the Bohr radius and the subscripts refer to quantum numbers $n$, $l$, and $m$, respectively. Sometimes these are written in terms of a dimensionless radius $\rho = r/a_0$. Sometimes $a_0$ is just explicitly set to $1$, though once this is done it's hard to figure out how to put dimensions back in, so I'll keep it there explicitly.
The $2P$ orbital is just the equally-weighted sum of these three parts:
$$ \psi_{2P} = \frac{1}{\sqrt{3}} \left(\psi_{2,1,0} + \psi_{2,1,1} + \psi_{2,1,-1}\right). $$
Here I put in the $\sqrt{3}$ in order to keep the probability normalized to $1$. If you are only interested in electron densities, you can take the square magnitude at this point, which I'll call $F$ to try to bring notations back in line. This turns out to be
$$ F(r, \theta, \phi) = \left\lvert \psi_{2P}(r, \theta, \phi) \right\rvert^2 = \frac{1}{96\pi a_0^5} r^2 \mathrm{e}^{-r/a_0} \left(1 - \sin^2(\theta) \cos(2\phi)\right), $$
as you are free to check.1
Now, to plot these, you are probably going to simply evaluate this function at Cartesian pixels.2 Therefore you want this as a function of $x$, $y$, and $z$. The usual definitions of spherical coordinates tell us
\begin{align}
r & = \sqrt{x^2+y^2+z^2}, \\
\sin^2(\theta) & = \frac{z^2}{x^2+y^2+z^2}, \\
\cos(2\phi) & = \frac{x^2-y^2}{x^2+y^2},
\end{align}
where the last identity is best seen by using the identity $\cos(2\phi) = \cos^2(\phi) - \sin^2(\phi)$.
At this point, simply plug the formulas back into the definition of $F$, setting one of $x$, $y$, or $z$ to $0$ at a time. For example, the $x{-}y$ slice corresponds to $z = 0$, so
$$ F_{x{-}y}(x, y) = \frac{x^2+y^2}{96\pi a_0^5} \mathrm{e}^{-\sqrt{x^2+y^2}/a_0} \left(1 - \frac{x^2-y^2}{(x^2+y^2)^2}\right). $$
That's all there is to the physics - the rest is just throwing this into the language of your choice.3
1 Students should check this, to get a better handle on efficiently spotting orthogonality relations when integrating well-chosen eigenfunctions of Hermitian operators. Only once you are old as me can you Mathematica your way out of this exercise ;)
2 That's certainly the standard behavior of imagesc(). My Matlab experience doesn't go much further - there are other stackexchanges for how to best code this.
3 In Matlab, a quick way would be something like
[X,Y] = meshgrid(-3:0.01:3);
F = (X.^2 + Y.^2) / (96*pi) * exp(-sqrt(X.^2+Y.^2)) * (1 - (X.^2-Y.^2)/(X.^2+Y.^2).^2);
imagesc(F)
