python 3.x – I have written this code of a Finite Element Method for solving a Boundary value problem. But when I run this code it does not give me any output

Here is the code for the Boundary value problem and when I write this code in Matlab with the Matlab syntax it gives me a graph. But when I run this python code it does not gives me any result. Can anyone figure out that what is wrong with this code or what I am missing here.

import numpy as np
    import matplotlib.pyplot as plt
    def BVP2():
        a=0        # interval [a,b]
        nNode = n+1  # number of nodes
        nElem = n     # number of elements
        x = np.zeros((nNode,1), dtype=int)
        # DOmain discretization
        for i in range(1, nNode):
            x = a + (i-1)*h 
        # Construction of finite elements
        for i in range(1, nElem):
            elem[i,:]=  [i,i+1]
            A = np.zeros((nNode,nNode))
            B = np.zeros((nNode,nNode))
            F = np.zeros((nNode,1))         
            u = np.zeros((nNode,1))
        for i in range(1, nElem):        
            node1 = elem(i,1)
            node2 = elem(i,2)
            Ae  = StiffMatrix1(x,h,node1,node2)
            Be  = StiffMatrix2(x,h,node1,node2)
            Fe  = Loadvector(x,h,node1,node2)
            A[node1:node2,node1:node2] = A[node1:node2,node1:node2] + Ae
            B[node1:node2,node1:node2] = B[node1:node2,node1:node2] + Be
            F[node1:node2] = F[node1:node2] + Fe
    # Implementation of BC's
            Ared  = A[2:nNode-1,2:nNode-1]
            Bred  = B[2:nNode-1,2:nNode-1]
            Fred = F[2:nNode-1]
            ured = u[2:nNode-1]         
            ured = [Ared + Bred]/Fred
            u = np.block([[0],[ured], [0]])
    # Exact Solution         
            uex = np.zeros((nNode,1))
            for i in range(1, nNode):
                uex = -(1/6)*x(i)**3 + (1/6)*x(i)             
            fig, ax = plt.subplots()
            ax.plot(x, u)
    def StiffMatrix1(x,h,node1,node2):
        x1 = x(node1)
        x2 = x(node2)
        A11 = (1/h)
        A12 = (-1/h)
        A21 = (-1/h)
        A22 = (1/h)
        Ae = np.block([[np.c_[A11, A12]], [np.c_[A21, A22]]])
        return Ae
    def StiffMatrix2(x,h,node1,node2):
        x1 = x(node1)
        x2 = x(node2)
        B11 = 0 #h/3
        B12 = 0 #h/6
        B21 = 0 #h/6
        B22 = 0 #h/3
        Be = np.block([[np.c_[B11, B12]], [np.c_[B21, B22]]])
        return Be
    def Loadvector(x,h,node1,node2):
        x1 = x(node1)
        x2 = x(node2)
        Fe1 = h/6*(x2+2*x1)
        Fe2 = h/6*(2*x2+x1)
        Fe = np.block([[Fe1],[Fe2]])
        return Fe

Your suggestions would be highly appreciated.

Read more here: Source link