FEM_Code := proc(  )
    uses ArrayTools, LinearAlgebra;
;
    #Finite element method code for solving bvp nonlinear ODEs%
    # u''+uu'-u=exp(2x), u(0)= 1, u(1)=e     %
    #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    clear("all");
    m_close("all");
    clc;
    n := 5;
    # NO of element
    nn := n +~ 1;
    # No of nodes
    lgth := 1;
    # Domain length
    he := MTM[mrdivide](lgth,n);
    # lenth of each elemnet
    x := Concatenate(2, Vector[row]([seq(0..lgth,he)]));
    # Data point for independant variable
    AC := 0.00005;
    # Accuracy
    F := rtable(proc(A)
    local j, n;
    if A :: rtable[row] then
       n := ArrayTools:-NumElems(A);
       if n = 1 then
           1..A(1), 1..A(1), fill = 0, subtype = Matrix;
       elif n = 2 then
           1..A(1), 1..A(2), fill = 0, subtype = Matrix;
       else
           seq(1..A(j), j = 1..n), fill = 0, subtype = Array;
       end if;
    else
       seq(1..args[i], i = 1..2), fill = 0, subtype = Array;
    end if;
end proc(nn, 1));
    # Initialization
    F(1) := evalhf(exp(0));
    F(nn) := evalhf(exp(1));
    # Boundary conditions
    #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    # Direct Iterative process to handle nonlinear problem
    #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    c := 1.0;
    count := 0;
    # Initializations for count for iterations
    tic;
    # Time start
    while ArrayTools[AllNonZero]((evalhf(zip(`>`,eval(c),eval(0))))) do
        F1 := assembly(F, n, he);
        c := 0.0;
        for Matlab_i from 1 to nn do
            if ArrayTools[AllNonZero]((evalhf(zip(`>`,eval(abs((F(Matlab_i) -~ F1(Matlab_i)))),eval(AC))))) then
                c := c +~ 1;
                break;
            end if;
        end do;
        F := copy(F1);
        count := count +~ 1;
    end do;
    print("Hence solution=:");
    #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    # Output for prinmary and secondary variables %%%
    #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    m_diff := abs((F -~ (HermitianTranspose(map[evalhf](exp,x)))));
    m_fprintf("No of element=%d\\", n);
    print("      x       FEM          Exact       Error");
    print(Concatenate(2, HermitianTranspose(x), F, HermitianTranspose(map[evalhf](exp,x)), m_diff));
    m_fprintf("No of iterations=%d\\", count);
    #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    #%% Ploting of primary variable %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    plots[display]([m_plot(x, F, "--rs"), plot("Linewidth", 2)]);
    xlabel("x");
    ylabel("u(x)");
    title("solution plot to given BVP");
    toc;
    #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    #%% Derivative of element matrix and Assembly%%%%%%%%%%%%%%%%%%%%
    #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    end proc;

assembly := proc( F, n, he )
    uses ArrayTools, LinearAlgebra;
    local F1;
    nn := n +~ 1;
    k := rtable(proc(A)
    local j, n;
    if A :: rtable[row] then
       n := ArrayTools:-NumElems(A);
       if n = 1 then
           1..A(1), 1..A(1), fill = 0, subtype = Matrix;
       elif n = 2 then
           1..A(1), 1..A(2), fill = 0, subtype = Matrix;
       else
           seq(1..A(j), j = 1..n), fill = 0, subtype = Array;
       end if;
    else
       seq(1..args[i], i = 1..2), fill = 0, subtype = Array;
    end if;
end proc(nn, nn));
    # Initialization of main Matrix
    R := rtable(proc(A)
    local j, n;
    if A :: rtable[row] then
       n := ArrayTools:-NumElems(A);
       if n = 1 then
           1..A(1), 1..A(1), fill = 0, subtype = Matrix;
       elif n = 2 then
           1..A(1), 1..A(2), fill = 0, subtype = Matrix;
       else
           seq(1..A(j), j = 1..n), fill = 0, subtype = Array;
       end if;
    else
       seq(1..args[i], i = 1..2), fill = 0, subtype = Array;
    end if;
end proc(nn, 1));
    # Initialization of RHS Matrix
    syms("x");
    # x as symbolic variable
    s := Concatenate(2, 1 -~ MTM[mrdivide](x,he), MTM[mrdivide](x,he));
    # linear shape function
    ds := m_diff(s, x);
    # Differentiations of shape function
    lmm := Vector();
    for Matlab_i from 1 to n do
        lmm := Concatenate(1, lmm, Concatenate(2, Matlab_i, Matlab_i +~ 1));
    end do;
    for Matlab_i from 1 to n do
        lm := lmm(Matlab_i, ..);
        #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        #%% Generation of Element Matrix k11 and RHS Matrix f1%%%%%%%%%%
        #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        k11 := -m_int(((HermitianTranspose(ds)) . ds), x, 0, he) +~ ((m_int((((HermitianTranspose(s)) . ds) . s(1)), x, 0, he) . F(lm(1))) +~ (m_int((((HermitianTranspose(s)) . ds) . s(2)), x, 0, he) . F(lm(2)))) -~ m_int(((HermitianTranspose(s)) . s), x, 0, he);
        f1 := m_int((map[evalhf](exp,(2 . (x +~ ((Matlab_i -~ 1) . he)))) . (HermitianTranspose(s))), x, 0, he);
        #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        #%% Assembly accroding to connectivity Matrix%%%%%%%%%%%%%%%%%%%%
        #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        k(lm, lm) := k(lm, lm) +~ k11;
        R(lm) := R(lm) +~ f1;
    end do;
    #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    #%% Imposing Boundary Conditions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    k(1, ..) := 0.0;
    k(nn, ..) := 0.0;
    k(1, 1) := 1.0;
    k(nn, nn) := 1.0;
    R(1, 1) := F(1);
    R(nn, 1) := F(nn);
    #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    #%% Solution of equations (F1) %%%%%%
    #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    d := MTM[mldivide](k,R);
    # better than using inverse k*R
    F1 := copy(d);
        return F1;
end proc;

with(ArrayTools):
with(LinearAlgebra):