import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

T = np.zeros([100,10])
for i in range(len(T)):
    T[i][-1] = 10

i = 0
while i < len(T)-1:
    for j in range(1,len(T[0])-1):
        T[i+1][j] = T[i][j] +(1./10)*(T[i][j-1] +T[i][j+1] -2*T[i][j])
    i += 1

plt.hold(True)
for i in range(len(T)):
    plt.plot(T[i])

np.linspace(0,np.pi,10)

T = np.zeros([100,20])
T[0] = np.sin(np.linspace(0,np.pi,20))**2

i = 0
while i < len(T)-1:
    for j in range(1,len(T[0])-1):
        T[i+1][j] = T[i][j] +(1./10)*(T[i][j-1] +T[i][j+1] -2*T[i][j])
    i += 1

plt.hold(True)
for i in range(len(T)):
    plt.plot(T[i])



T = np.zeros([100,20])
T[0][1:-1] = 100.
T

i = 0
while i < len(T)-1:
    for j in range(1,len(T[0])-1):
        T[i+1][j] = T[i][j] +(1./10)*(T[i][j-1] +T[i][j+1] -2*T[i][j])
    i += 1

plt.hold(True)
for i in range(len(T)):
    plt.plot(T[i])



from ipywidgets import StaticInteract, RangeWidget, RadioWidget

def update(T,t):
    for i in range(0,t):
        for j in range(1,len(T[0])-1):
            T[i+1][j] = T[i][j] +(1./10)*(T[i][j-1] +T[i][j+1] -2*T[i][j])
    return T[t]
    
def plot(t):
    T = np.zeros([10,10])
    for i in range(len(T)):
        T[i][-1] = 10
    fig, ax = plt.subplots(figsize=(4, 3),
                           subplot_kw={'axisbg':'#EEEEEE',
                                       'axisbelow':True})
    ax.grid(color='w', linewidth=2, linestyle='solid')
    array = update(T,t)
    ax.plot(array, lw=5, alpha=0.4, label = t)
    ax.legend(loc='upper right')
    return fig

StaticInteract(plot, t=RangeWidget(1., 4., 1.))

for t in range(len(T)):
    plt.hold(True)
    array = update(T,t)
    plt.plot(array)

T = np.zeros([10,10])
T[0][1:-1] = 100

op = np.zeros([8,8])
a = 1./2
for i in range(len(op)):
    for j in range(len(op[0])):
        if i==j:
            op[i][j] = 1+2*a
        elif abs(i-j) == 1:
            op[i][j] = -a

op

opinv = np.linalg.inv(op)

len(opinv)

for i in range(1,len(T)):
    T[i][1:-1] = np.dot(opinv,T[i-1][1:-1])

plt.hold(True)
for i in range(len(T)):
    plt.plot(T[i])

len(T[0][1:-1])

I = np.identity(8)
B = np.zeros([8,8])
a = 1./2
for i in range(len(B)):
    for j in range(len(B[0])):
        if i==j:
            B[i][j] = 2
        elif abs(i-j) == 1:
            B[i][j] = -1

op = np.dot(np.linalg.inv(2*I+a*B),2*I-a*B)

for i in range(1,len(T)):
    T[i][1:-1] = np.dot(opinv,T[i-1][1:-1])

plt.hold(True)
for i in range(len(T)):
    plt.plot(T[i])

from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = plt.axes(projection='3d')
#ax.scatter(x,y,T)
ax.plot_surface(x,y,T)

x = np.arange(len(T))
y = np.arange(len(T))
x,y = np.meshgrid(x,y)