% fem.m
% Solves a transient diffusion problem using the
% finite element method with linear basis functions.
% needs lu_tri and fas_tri
global nt delx2 afd bfd cfd

% set up the parameters
nt = 17
delx = 1/(nt-1)
delx2 = delx*delx
x(1) = 0.;
for i=2:nt-1
c0(i-1) = 0.;
x(i) = x(i-1)+delx;
end
x(nt) = x(nt-1)+delx

% do the LU decomposition of the CC matrix
for i=1:nt-2
afd(i) = 1/6;
bfd(i) = 2/3;
cfd(i) = 1/6;
end
lu_tri(nt-2)

tbreak = [0 .001 .01 .1 1];
for kkk=1:4
tspan = [tbreak(kkk) tbreak(kkk+1)];
[t cc] = ode45('fem_rhs',tspan,c0);
nn=size(t);
c(kkk,:) = cc(nn(1),:);
c0(:) = c(kkk,:);
end

% add in the first and last point
for kkk=1:4
d(kkk,1) = 1;
for i=2:nt-1
d(kkk,i)=c(kkk,i-1);
end
d(kkk,nt) = 0;
end

plot(x,d(4,:),'g-*',x,d(3,:),'m-*',x,d(2,:),'b-*',x,d(1,:),'r-*')
xlabel('x')
ylabel('c')
legend('t = 1.0','t = 0.1','t = 0.01','t = 0.001')
title('Diffusion Problem Solved using Galerkin Finite Element Method with Linear Basis Functions')

******************************************************************
% fem_rhs.m
function ydot=fem_rhs(time,y)
global nt delx2 afd bfd cfd dfd

% add in the boundary conditions
c(1) = 1;
c(nt) = 0;

% transfer the y to the c
for i=1:nt-2
c(i+1) = y(i);
end

% calculate the right-hand side
for i=1:nt-2
dfd(i) = c(i) - 2*c(i+1) + c(i+2);
end
dfd = dfd/delx2;

% call the fore and aft sweep
fas_tri(nt-2)

% put result in the ydot
ydot = dfd';

******************************************************************
function lu_tri(n)
global afd bfd cfd

% this subroutine does an lu decomposition of
% a tridiagonal system of
% equations of the type that often occur with the
% finite difference method.
% After calling lu_tri one must call fas_tri to
% process the right-hand side.
% input
% a(n), b(n), c(n) in the following equation:
% a(i)*y(i-1) + b(i)*y(i) + c(i)*y(i+1) = d(i)
% n - the size of the system being solved.

% output
% afd(n), bfd(n), cfd(n) this is the lu decomposition.
% the original afd, bfd, cfd are destroyed.

% lower decomposition tridiagonal matrix

for l=2:n
s = afd(l)/bfd(l-1);
bfd(l) = bfd(l)-s*cfd(l-1);
afd(l) = s;
end

******************************************************************

function fas_tri(n)
global afd bfd cfd dfd

% this subroutine does the fore and aft
% sweep to solve a tridiagonal system of
% equations of the type that often occur with the
% finite difference method.
% One must call lu_tri one must call first.
% input
% afd(n), bfd(n), cfd(n) from the lu decomposition
% solving the problem:
% a(i)*y(i-1) + b(i)*y(i) + c(i)*y(i+1) = d(i)
% n - the size of the system being solved.
% output
% dfd(n) is the solution

% forward sweep tridiagonal matrix
for l=2:n
dfd(l) = dfd(l)-afd(l)*dfd(l-1);
end
% back substitution
dfd(n) = dfd(n)/bfd(n);
for l=2:n
k = n-l+1;
dfd(k) = (dfd(k)-cfd(k)*dfd(k+1))/bfd(k);
end