The Wrong Code Will often Provide Beautiful Result, :)

 It means to compute 2d diffusion equation just like previous post in polar/cylindrical coordinate, and all went to wrong direction, :)

 Still trying to understand matplotlib mplot3d behavior

import scipy as sp

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.pyplot as plt

import mpl_toolkits.mplot3d.axes3d as p3
import matplotlib.animation as animation

#dr  = .1
#dp  = .1
#nr      = int(1/dr)
#np      = int(2*sp.pi/dp)

nr  = 10
np  = 10

dr  = 1./nr
dp  = 2*sp.pi/np

a   = .5

tmax    = 100
t       = 0.


dr2     = dr**2
dp2     = dp**2

dt      = dr2 * dp2 / (2 * a * (dr2 + dp2) )
dt      /=10.
print dt



ut      = sp.zeros([nr,np])
u0      = sp.zeros([nr,np])

ur      = sp.zeros([nr,np])
ur2     = sp.zeros([nr,np])

r       = sp.arange(0.,1.,dr)
p       = sp.arange(0.,2*sp.pi,dp)

#initial

for i in range(nr):
    for j in range(np):
        if ( (i>(2*nr/5.)) & (i<(3.*nr/3.)) ):
            u0[i,j] = 1.
#print u0

def hitung_ut(ut,u0):
    for i in sp.arange (len(r)):
        if r[i]!= 0.:
            ur[i,:]     = u0[i,:]/r[i]
            ur2[i,:]     = u0[i,:]/(r[i]**2)
    ut[1:-1, 1:-1]  = u0[1:-1, 1:-1] + a*dt*(
            (ur[1:-1, 1:-1] - ur[:-2, 1:-1])/dr+
            (u0[2:, 1:-1] - 2*u0[1:-1, 1:-1] + u0[:-2,1:-1])/dr2+
            (ur2[1:-1, 2:] - 2*ur2[1:-1, 1:-1] + ur2[1:-1, :-2])/dp2)


#hitung_ut(ut,u0)
#print ut


def data_gen(framenumber, Z ,surf):
    global ut
    global u0
    hitung_ut(ut,u0)
    u0[:] = ut[:]
    Z = u0
    
    ax.clear()
    plotset()
    surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False, alpha=0.7)
    return surf,


fig = plt.figure()
#ax = fig.gca(projection='3d')
ax = fig.add_subplot(111, projection='3d')

R = sp.arange(0,1,dr)
P = sp.arange(0,2*sp.pi,dp)
R,P = sp.meshgrid(R,P)

X,Y = R*sp.cos(P),R*sp.sin(P) 

Z = u0
print len(R), len(P)



def plotset():
    ax.set_xlim3d(-1., 1.)
    ax.set_ylim3d(-1., 1.)
    ax.set_zlim3d(-1.,1.)
    ax.set_autoscalez_on(False)
    ax.zaxis.set_major_locator(LinearLocator(10))
    ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
    cset = ax.contour(X, Y, Z, zdir='x', offset=0. , cmap=cm.coolwarm)
    cset = ax.contour(X, Y, Z, zdir='y', offset=1. , cmap=cm.coolwarm)
    cset = ax.contour(X, Y, Z, zdir='z', offset=-1., cmap=cm.coolwarm)

plotset()
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False, alpha=0.7)

fig.colorbar(surf, shrink=0.5, aspect=5)

ani = animation.FuncAnimation(fig, data_gen, fargs=(Z, surf),frames=500, interval=30, blit=False)
#ani.save('2dDiffusionf500b512.mp4', bitrate=512)

plt.show()    

.