Example: solving ODE using FEM method

clc;
clear;
%  find the solution of du/dt -1 =0 Using finite element method
% the zone of solution is [0,1]
E = 5; % the number of elements
N = E+1; % the number of nodes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%creat the gids
h = 1.0/E;
element = zeros(5,2); % the elements data
for i = 1:5
    element(i,1) = i;
end
for i = 1:5
    element(i,2) =i+1;
end

node = 0:h:1.0; % the nodes data
boundary = [0;0];% the boundary value
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create element matrix
ke = [-0.5 0.5;-0.5 0.5;];
be = h*[0.5;0.5];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%create main matrix
k = zeros(N,N);
b = zeros(N,1);
%%%%
for i = 1:E
    for m = 1:2
        for n = 1:2
            k(element(i,m),element(i,n)) = k(element(i,m),element(i,n)) + ke(m,n);
        end
    end
end
for i  = 1:E
    for m =1:2
        b(element(i,m),1)= b(element(i,m),1) + be(m,1);
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% consider the boundary value u(0) = 0  in node = 1
b = b- boundary(1)*k(:,1);
for i = 1:N
    k(1,i) = 0;
    k(i,1) = 0;
end
k(1,1) = 1.0;
u = k\b;
u(1) = boundary(1);
u;