Okay so i was able to accomplish what i wanted to in Cython, thank god.
In the process I found Cython much more powerful of a tool !

So the result is a 2D wave PDE solved using a finite difference. The
main computation has been exported to a Cythonized function. The function
takes in three 2D np.ndarrays and returns one back. Much easier to work
with numpy and other data types as well.

Here is the Cythonized function:

```
from numpy import *
cimport numpy as np
def cwave_prop( np.ndarray[double,ndim=2] u_prev, np.ndarray[double,ndim=2]
u, np.ndarray[double,ndim=2] u_next):
cdef double t = 0
cdef double t_old = 0
cdef double t_end = 100
cdef int i,j
cdef double c = 1
cdef double Lx = 100
cdef double Ly = 100
cdef double dx = 1
cdef double dy = 1
cdef double dt = (1/(c))*(1/(sqrt(1/dx**2 + 1/dy**2)))
while t<t_end:
t_old=t; t +=dt
#wave steps through time and space
for i in xrange(1,99):
for j in xrange(1,99):
u_next[i,j] = - u_prev[i,j] + 2*u[i,j] +
((c*dt/dx)**2)*(u[i-1,j] - 2*u[i,j] + u[i+1,j]) +
((c*dt/dx)**2)*(u[i,j-1] - 2*u[i,j] + u[i,j+1])
#set boundary conditions of grid to 0
for j in xrange(100): u_next[0,j] = 0
for i in xrange(100): u_next[i,0] = 0
for j in xrange(100): u_next[Lx-1,j] = 0
for i in xrange(100): u_next[i,Ly-1] = 0
#set prev time step equal to current one
for i in xrange(100):
for j in xrange(100):
u_prev[i,j] = u[i,j];
u[i,j] = u_next[i,j];
print u_next
```

And here is my python script that calls it and also returns it back then
plots the result.
Any suggestions on writing better code are welcome...

```
'''George Lees Jr.
2D Wave pde '''
from numpy import *
import numpy as np
import matplotlib.pyplot as plt
import cwave2
np.set_printoptions(threshold=np.nan)
#declare variables
#need 3 arrays u_prev is for previous time step due to time derivative
Lx=Ly = (100) #Length of x and y dims of grid
dx=dy = 1 #derivative of x and y respectively
x=y = np.array(xrange(Lx)) #linspace to set the initial
condition of wave
u_prev=np.ndarray(shape=(Lx,Ly), dtype=np.double) #u_prev 2D grid for
previous time step needed bc of time derivative
u=np.ndarray(shape=(Lx,Ly), dtype=np.double) #u 2D grid
u_next=np.ndarray(shape=(Lx,Ly), dtype=np.double) #u_next for advancing
the time step #also these are all numpy ndarrays
c = 1 #setting constant velocity of the wave
dt = (1/float(c))*(1/sqrt(1/dx**2 + 1/dy**2)) #we have to set dt
specifically to this or numerical approx will fail!
print dt
#set Initial Conditions and Boundary Points
#I(x) is initial shape of the wave
def I(x,y): return exp(-(x-Lx/2.0)**2/2.0 -(y-Ly/2.0)**2/2.0)
#set up initial wave shape
for i in xrange(100):
for j in xrange(100):
u[i,j] = I(x[i],y[j])
#set up previous time step array
for i in xrange(1,99):
for j in xrange(1,99):
u_prev[i,j] = u[i,j] + 0.5*((c*dt/dx)**2)*(u[i-1,j] - 2*u[i,j]
+ u[i+1,j]) +
0.5*((c*dt/dy)**2)*(u[i,j-1] - 2*u[i,j] + u[i,j+1])
#set boundary conditions to 0
for j in xrange(100): u_prev[0,j] = 0
for i in xrange(100): u_prev[i,0] = 0
for j in xrange(100): u_prev[Lx-1,j] = 0
for i in xrange(100): u_prev[i,Ly-1] = 0
#call C function from Python
cwave2.cwave_prop( u_prev , u , u_next )
#returned u (2D np.ndarray)
from tempfile import TemporaryFile
outfile = TemporaryFile()
np.save(outfile,u)
fig = plt.figure()
plt.imshow(u,cmap=plt.cm.ocean)
plt.colorbar()
plt.show()
```