import numpy as np
import matplotlib.pyplot as plt

def fea_solve(Kpp, Kpf, Kfp, Kff, xp, Ff, M):
    xf = my_triangular_solve(M, Ff - np.matmul(Kfp, xp))
    return xf


M = my_lu(Kff)
xf_1 = fea_solve(Kpp, Kpf, Kfp, Kff, xp, Ff[:, 0], M)
xf_2 = fea_solve(Kpp, Kpf, Kfp, Kff, xp, Ff[:, 1], M)

image_xf_1 = plot_truss(xf_1)
image_xf_2 = plot_truss(xf_2)