%Finite element method code for solving bvp nonlinear ODEs%
% u''+uu'-u=exp(2x), u(0)= 1, u(1)=e     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function FEM_Code()
clear all; close all; clc
n=5;                     % NO of element
nn=n+1;                  % No of nodes
lgth=1;                  % Domain length
he=lgth/n;               % lenth of each elemnet
x=[0:he:lgth];           % Data point for independant variable
AC=0.00005;              % Accuracy
F=zeros(nn,1);           % Initialization
F(1)=exp(0); F(nn)=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 (c>0)
        [F1]=assembly(F,n,he);
          c=0.0;
          for i=1:nn
            if (abs(F(i)-F1(i))>AC)
                 c=c+1;
                 break;
            end
          end 
         F=F1;
         count=count+1;
end
  disp('Hence solution=:');
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Output for prinmary and secondary variables %%%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  diff=abs(F-exp(x)');
  fprintf('No of element=%d\n',n)
  disp('      x       FEM          Exact       Error')
  disp([x',F,exp(x)',diff])
  fprintf('No of iterations=%d\n',count)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%% Ploting of primary variable %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  plot(x,F,'--rs','Linewidth',2)
  xlabel('x')
  ylabel('u(x)')
  title('solution plot to given BVP')
  toc                                % given totlal time
  
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Derivative of element matrix and Assembly%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [F1]=assembly(F,n,he)
nn=n+1;
k = zeros(nn,nn);               % Initialization of main Matrix
R = zeros(nn,1);                % Initialization of RHS Matrix
syms x                          % x as symbolic variable
s=[1-x/he,x/he];                % linear shape function
ds=diff(s,x);                   % Differentiations of shape function
lmm =[];
for i=1:n
    lmm=[lmm;[i,i+1]];          % connectvity Matrix
end
for i=1:n
    lm=lmm(i,:);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% Generation of Element Matrix k11 and RHS Matrix f1%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    k11=-int(ds'*ds,x,0,he)+(int(s'*ds*s(1),x,0,he)*F(lm(1))...
        +int(s'*ds*s(2),x,0,he)*F(lm(2)))-int(s'*s,x,0,he);
    f1 = int(exp(2*(x+(i-1)*he))*s',x,0,he);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% Assembly accroding to connectivity Matrix%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    k(lm,lm) = k(lm,lm) + k11;
    R(lm) = R(lm) + f1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 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 = k\R;           % better than using inverse k*R
F1 = d;
end

    

                                
 


  
        
      






;