For a project I am working on (in hyperbolic PDEs) I would like to get some rough handle on the behavior by looking at some numerics. I am, however, not a very good programmer. Can you recommend some resources for learning how to effectively code finite difference schemes in Scientific Python (other languages with small learning curve also welcome)?To give you an idea of the audience (me) for this recommendation:I am a pure mathematician by training, and am somewhat familiar with the theoretical aspects of finite difference schemesWhat I need help with is how to make the computer compute what I want it to compute, especially in a way that I don't duplicate too much of the effort already put in by others (so as to not re-invent the wheel when a package is already available). (Another thing I would like to avoid is to stupidly code something by hand when there are established data structures fitting the purpose.)I have had some coding experience; but I have had none in Python (hence I don't mind if there are good resources for learning a different language [say, Octave for example]). Books, documentation would both be useful, as would collections of example code. 
Recommendation for Finite Difference Method in Scientific Python
python;finite difference;reference request;hyperbolic pde
Here is a 97-line example of solving a simple multivariate PDE using finite difference methods, contributed by Prof. David Ketcheson, from the py4sci repository I maintain.  For more complicated problems where you need to handle shocks or conservation in a finite-volume discretization, I recommend looking at pyclaw, a software package that I help develop.Pattern formation code    Solves the pair of PDEs:       u_t = D_1 \nabla^2 u + f(u,v)       v_t = D_2 \nabla^2 v + g(u,v)import matplotlibmatplotlib.use('TkAgg')import numpy as npimport matplotlib.pyplot as pltfrom scipy.sparse import spdiags,linalg,eyefrom time import sleep#Parameter valuesDu=0.500; Dv=1;delta=0.0045; tau1=0.02; tau2=0.2; alpha=0.899; beta=-0.91; gamma=-alpha;#delta=0.0045; tau1=0.02; tau2=0.2; alpha=1.9; beta=-0.91; gamma=-alpha;#delta=0.0045; tau1=2.02; tau2=0.; alpha=2.0; beta=-0.91; gamma=-alpha;#delta=0.0021; tau1=3.5; tau2=0; alpha=0.899; beta=-0.91; gamma=-alpha;#delta=0.0045; tau1=0.02; tau2=0.2; alpha=1.9; beta=-0.85; gamma=-alpha;#delta=0.0001; tau1=0.02; tau2=0.2; alpha=0.899; beta=-0.91; gamma=-alpha;#delta=0.0005; tau1=2.02; tau2=0.; alpha=2.0; beta=-0.91; gamma=-alpha; nx=150;#Define the reaction functionsdef f(u,v):    return alpha*u*(1-tau1*v**2) + v*(1-tau2*u);def g(u,v):    return beta*v*(1+alpha*tau1/beta*u*v) + u*(gamma+tau2*v);def five_pt_laplacian(m,a,b):    Construct a matrix that applies the 5-point laplacian discretization    e=np.ones(m**2)    e2=([0]+[1]*(m-1))*m    h=(b-a)/(m+1)    A=np.diag(-4*e,0)+np.diag(e2[1:],-1)+np.diag(e2[1:],1)+np.diag(e[m:],m)+np.diag(e[m:],-m)    A/=h**2    return Adef five_pt_laplacian_sparse(m,a,b):    Construct a sparse matrix that applies the 5-point laplacian discretization    e=np.ones(m**2)    e2=([1]*(m-1)+[0])*m    e3=([0]+[1]*(m-1))*m    h=(b-a)/(m+1)    A=spdiags([-4*e,e2,e3,e,e],[0,-1,1,-m,m],m**2,m**2)    A/=h**2    return A# Set up the grida=-1.; b=1.m=100; h=(b-a)/m; x = np.linspace(-1,1,m)y = np.linspace(-1,1,m)Y,X = np.meshgrid(y,x)# Initial datau=np.random.randn(m,m)/2.;v=np.random.randn(m,m)/2.;plt.hold(False)plt.pcolormesh(x,y,u)plt.colorbar; plt.axis('image'); plt.draw()u=u.reshape(-1)v=v.reshape(-1)A=five_pt_laplacian_sparse(m,-1.,1.);II=eye(m*m,m*m)t=0.dt=h/delta/5.;plt.ion()#Now step forward in timefor k in range(120):    #Simple (1st-order) operator splitting:    u = linalg.spsolve(II-dt*delta*Du*A,u)    v = linalg.spsolve(II-dt*delta*Dv*A,v)    unew=u+dt*f(u,v);    v   =v+dt*g(u,v);    u=unew;    t=t+dt;    #Plot every 3rd frame    if k/3==float(k)/3:        U=u.reshape((m,m))        plt.pcolormesh(x,y,U)        plt.colorbar        plt.axis('image')        plt.title(str(t))        plt.draw()plt.ioff()