clear all;

% Integral Limits
Tmin=0;
Tmax=1;
% Number of points in temperature
K=100;
% Increment of temperature
dT=(Tmax-Tmin)/K;

T(1)=Tmin;
c(1)=0;  %heat capacity at T=0

% number of points used in Simpson rule.
N=100;

for j=2:K+1
    
    % Temperature
    T(j)=dT*(j-1)
    
    % Step size
    h=(T(j)-Tmin)/N;
    
    % value at the lower limit
    x(1)=Tmin;
    f(1)=0;
    
    % generate the function values
    for i=1:N
       x(i+1)=Tmin+i*h;
       f(i+1)=exp(-1/x(i))/(1+exp(-1/x(i)))^2;
    end
    
    % Simpson rule
    c(j)=sum(f(1:2:N-1)+4*f(2:2:N)+f(3:2:N+1))*h/3;

end
    

p=plot(T,c);
xlabel('T');
ylabel('Correction Term');
set(p(1),'Color','blue','Linewidth',2);
