Plotting a conditional function using matplotlib in Python
<p>I'm trying to create a line or scatter plot of this algorithm and it gives me the error </p>

<blockquote>
  <p>Traceback (most recent call last):
  File "/Users/itstest/Documents/workspace/Practice/src/PlutoModel.py", line 73, in module 
  plt.plot(xr, P(xr))
  File "/Users/itstest/Documents/workspace/Practice/src/PlutoModel.py", line 55, in P
    if x > r:
  ValueError: "The truth value of an array with more than one element is ambiguous. Use <code>a.any()</code> or <code>a.all()</code>."</p>
</blockquote>

<p>I have looked up possible solutions to this error but I don't think any apply to me.</p>

<pre><code>import numpy as np
import scipy.integrate as integ
import matplotlib.pyplot as plt

rho = 1860
rhos = 250 #Assuming Nitrogen Ice
rhom = 1000 #Assuming Water
rhoc = 3500 #Assuming a mix of Olivine and Pyroxene

def rhof(x):
    if x &gt; r:
        return "Point is outside of the planet"
    elif x &lt; r and x &gt; rm:
        return rhos
    elif x &lt; rm and x &gt; rc:
        return rhom
    else:
        return rhoc

r = 1.187e6
rc = 8.5e5 #Hypothesized
rm = 9.37e5 #Estimated based on crustal thickness of 250 km

Ts = 44
B = 0.15
G = 6.67e-11

m = 1.303e22
mc = (4*np.pi*rhoc*rc**3)/3
mm = (4*np.pi*rhom*((rm**3) - (rc**3)))/3
ms = (4*np.pi*rhos*((r**3) - (rm**3)))/3

Ic = 0.4*mc*rc**2
Im = 0.4*mm*rm**2
Is = 0.4*ms*r**2
Itot = Is + Im + Ic

def gi(x):
    if x == r:
        return G*m/r**2
    elif x &gt; r:
        return "Point is outside of the planet"
    elif x &gt; rc and x &lt; rm:
        return (G*mc/rc**2) + (G*mm/((x**2) - (rc**2)))
    elif x &gt; rm and x &lt; r:
        return (G*mc/rc**2) + (G*mm/((rm**2) - (rc**2))) + (G*ms/((x**2) - (rm**2)))
    else:
        return G*((3*rhoc)/4*np.pi*x**3)/x**2   

def Psmb(z):
    return rhos*G*(4.0/3.0)*np.pi*(1/z**2)*(rhom*(rm**3) + rhos*(z - rm**3))
def Pmcb(z):
    return rhom*G*(4.0/3.0)*np.pi*(1/z**2)*(rhoc*(rc**3) + rhom*(z - rc**3))
def P(x):
    if x &gt; r:
        return "The point is outside of the planet"
    elif x == r:
        return 1
    elif x &gt; rm and x &lt; r:
        return (integ.quad(1000*gi(x), x, r))[0]
    elif x == rm:
        return (integ.quad(Psmb, x, r))[0]
    elif x &gt; rc and x &lt; rm:
        return (integ.quad(1000*gi(x), x, rm) + P(rm))[0]
    elif x == rc:
        return (integ.quad(Pmcb, x, rm) + P(rm))[0]
    elif x &lt; rc and x != 0:
        return (integ.quad(1000*gi(x), x, rc) + P(rc))[0]
    else:
        return ((2.0/3.0)*np.pi*G*(rhoc**2)*r**2)

xr = np.linspace(0, 1187000, 1000)
plt.plot(xr, P(xr))

print("Mass = " + str(m) + " kg")
print("Radius = " + str(r) + " m")
print("Density = " + str(rho) + " kg/m^3")
print("Moment of Inertia = " + str(Itot) + " kgm^2")
print("Mean Surface Temperature = " + str(Ts) + " K")
print("Magnetic Field = " + str(B) + " nT")
print("Surface Gravity = " + str(gi(r)) + " m/sec^2")
print("Pressure at Surface = " + str(P(r)) + " Pa")
print("Pressure at Crust-Mantle Boundary = " + str(P(rm)) + " Pa")
print("Pressure at Mantle-Core Boundary = " + str(P(rc)) + " Pa")
print("Pressure at the Center = " + str(P(0)) + " Pa")
</code></pre>

<p>Is there a way to plot this function without separating each condition?</p>