Simple 1d implicit Euler Finite Difference Heat flow calculation

1d-heat-df.py

import numpy as np

n = 10
dt = 1
alpha = 1e-6 * 3.15e7
dx = 10
c = dt * alpha / (dx)**2
z = np.arange(0,n*dx,dx)
Q = 0 # change this
k = 1 # change this too

# Define FD scheme
B = np.diag([1+2*c] * n) + np.diag([-c] * (n-1),k=1)\
    + np.diag([-c] * (n-1),k=-1)

un = np.zeros((n,1))
unm1 = np.zeros((n,1))
unm1[int(n/2)] = 1

''' 
Set up time loop here
'''

un = np.linalg.solve(B,unm1)

# Enforce Dirichlet boundary conditions
un[0] = 1
# un[-1] = 0

# FD coefficients
c1 = -4
c2 = 1
c3 = -3

un[-1] = (c1*un[-2] + c2*un[-3] - Q/k * 2*dx)/c3

import matplotlib.pyplot as plt
plt.plot(z,unm1)
plt.plot(z,un)