Some answers are starting to get there... But the points are being connected into a line on the plot. How do we just plot the points?

import matplotlib.pyplot as plt
import numpy as np

def f(x):
 if(x == 2): return 0
 else: return 1

x = np.arange(0., 5., 0.2)

y = []
for i in range(len(x)):
   y.append(f(x[i]))

print x
print y

plt.plot(x,y,c='red', ls='', ms=5, marker='.')
ax = plt.gca()
ax.set_ylim([-1, 2])

plt.show()

#!/usr/bin/env python
import numpy as np
import matplotlib.pyplot as plt 

def fun (n, x): 
    if n <= x <= n + 1:
        return float(x) - n 
    elif n + 1 <= x <= n + 2:
        return 2.0 - x + n 
    return 0.0 

vfun = np.vectorize(fun)

x = np.linspace(0, 10, 1000)    
y = vfun(3, x)

plt.plot(x, y, '-')
plt.show()

There are quite a few issues here:

• None of the keyword arguments (Constant, Variable, unit, Real) that you are passing to Symbol are things that are recognized by SymPy. The only one that is close is real, which should be lowercase (like Symbol('x', real=True)). The rest do nothing. If you want units, you should use the SymPy units module in sympy.physics.units. There is no need to specify if a symbol is constant or variable.

• You have redefined Solar_Radius and gamma as numbers. That means that the Symbol definitions for those variables are pointless.

• If you are using Python 2, make sure to include from __future__ import division at the top of the file, or else things like 1/2 and 5/3 will be truncated with integer division (this isn't an issue in Python 3).

• range(0, x/c) doesn't make sense. range creates a list of numbers, like range(0, 3) -> [0, 1, 2]. But x/c is not a number, it's a symbolic expression.

• Additionally, xi[x, t] = ... doesn't make sense. xi is a Symbol, which doesn't allow indexing and certainly doesn't allow assignment.

• Don't mix numeric (math, mpmath, numpy, scipy) functions with SymPy functions. They won't work with symbolic expressions. You should use only SymPy functions. If you create a symbolic expression and want to convert it to a numeric one (e.g., for plotting), use lambdify.

What you want here is Piecewise. The syntax is Piecewise((expr, cond), (expr, cond), ..., (expr, True)), where expr is an expression that is used when cond is true ((expr, True) is the "otherwise" condition).

For your example, I believe you want

expr = Piecewise((0, t < x/c), (sy.exp(gamma*g*x/(2*c**2))*sy.besselj(0, (gamma*g/(2*c)*sy.sqrt(t**2 - (x/c)**2)))/2, t >= x/c))

If you want to turn this into a numeric function in x and t, use

xi = lambdify((x, t), expr)