Miscellany

import numpy as np
#r_i=7
#r_k=5
maxtime=3.0
l=1.0
#r_k=20
#r_i=40
eta=1.0
rho=1.0
omega=0.1
#ht=maxtime/r_k
#hx=l/r_i
hx=0.1
ht=0.1
K=maxtime/ht
I=l/hx
r_k=int(K)
r_i=int(I)
a=eta/rho
t=3

u=range(r_i)
for i in range(r_i):
	u[i]=range(r_k)
#for a in range(i):
#	for b in range(k):
#		arr[a][b]=0
def phi(x,t):
	if 0<=x<=l/3:
		return (1/rho)*(3*x/l)*np.sin(omega*t)
	else:
		return (1/rho)*(-1.5*(x/l-1))*np.sin(omega*t)
		
for i in range(r_i):
	u[i][0]=0
	u[i][1]=0
for i in np.arange(1, r_i-1):
	for k in np.arange(1, r_k-1):
		u[i][k+1]=a**2*ht**2**(u[i+1][k]-2*u[i][k]+u[i-1][k])/hx**2+phi(i*hx, t)*ht**2+2*u[i][k]-u[i][k-1]
for i in range(r_i):
	for k in range(r_k):
		print u[i][k],
	print