In [1]:
from IPython.display import display
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

In [2]:
k1 = 10
k2 = 1
a = 0
b = 0.5
c = 1
ua = 0
uc = 1
f = 1

In [3]:
h = 0.1
n = int((c - a) / h)
m = int((b-a)/c * n)-1

In [4]:
def mat_f(i, j):
    if i == m:
        if i == j:
            return k1 + k2
        if j == i-1:
            return -k1
        if j == i+1:
            return -k2
    if i < m:
        cons = k1
    elif i > m:
        cons = k2
    
        
    if i == j:
        return 2*cons
    elif j == i-1 or j == i+1:
        return -cons
    else:
        return 0

In [5]:
M = np.fromfunction(np.vectorize(mat_f), (n-1, n-1))
M *= 1/h

In [6]:
y = h * f * np.ones(n-1)
y[0] += k1 * 1/h * ua
y[n-2] += k2 * 1/h * uc

In [7]:
x = np.dot(np.linalg.inv(M),y)
u = np.hstack((np.array([ua]), x, np.array([uc])))

In [8]:
plt.plot(u)
plt.show()



In [9]:
koef1 = k1 * np.ones(m+1)
koef2 = k2 * np.ones(n-m-1)
koeff = np.hstack([koef1, koef2])
udot = np.diff(u)/h
plt.plot(udot*koeff)
plt.show()



In [ ]: