clear, clc, close all, clf
 fluxleft = 10^12;  Ea = 1.25;  D = 1/(3*Ea); L= sqrt(D/Ea);  h = 0.01; x = [0:.01:2.5 + 2*D];
 %I evaluated it from 0 to 2.5 + 2*D because the flux for the anaytical
 %solution must vanish at the extrapolation distance. In order to plot
 %relative error, the numerical solution must also be evaluated at these
 %bounds and is as to why the input parameters for the right boundry
 %condition was altered. The overall effect is negligible for when you 
 %are to evaluate flux from 0 to 2.5. (Atleast I hope so.) If this is a
 %problem an easy fix is to simply take the bound from 0 to 2.5.
 
  M =eye(length(x))*(2*D/h^2 + Ea);
for i = 1:(length(x)-2)
    M(i,i+1) = -D/h^2;
    M(i + 1,i) = -D/h^2;
end
M(1,1) = D/h^2 + 2*D/h^2 + Ea;
M(1,2) = -D/h^2;
M(end,end) = D/h^2 + 1/2.1312/h + Ea;
M(end,end-1) = -D/h^2;
 M(end-1,end) = -D/h^2;

a = zeros(length(x)); S = a(:,1);
S(1) = 2*D*fluxleft/h^2;

FLUX = M\S;
plot(x,FLUX, 'g'); grid on; grid minor; xlabel('\bf x'); ylabel('\it \bf {\0}(x)'); title('Attenuation of {\0}(x)'); hold on
x2 = [-.5:.01:0];
plot(x2, fluxleft*ones(1,length(x2)))
xline(0); xline(2.5), yline(0); xline(2.5+2*D); xlim([-.5, 2.5+2*D+0.2]);



y =   fluxleft/(sinh(-1/L*(2.5 + 2*D)))*sinh(1/L*(x - 2.5 - 2*D)); %Ideal Solution
y2 = 10^12*exp(-1/L*x);
plot(x,y, 'b')
plot(x,y2, 'k')
ylim([0,  fluxleft/10 + fluxleft])
legend('Numerical','Initial Flux', 'Analytical (sinh)', 'Analytical (exp)')


for i = 1:length(x)
relerr(i) = abs(1 - FLUX(i)/y(i));
relerr2(i) = abs(1 - FLUX(i)/y2(i));
end
figure(2)
 plot(x,relerr); hold on
 plot(x,relerr2, 'k'); xlim([0 2.5])
 legend( 'Sinh', 'Exp')
 grid on; grid minor; xlabel('\bf x'); ylabel('\bf Relative Error (decimal)'); title('Relative Error as a Function of x'); hold on
 FLUX(50)
 FLUX(100)