Heat equation analytical solution

heat.m

%% initialization

a=[1 10 100];
N=100;
h=2*pi/N;
T=1;

%time discretization is given by requirements for numerical stability
n=floor(log(h^2./(2*a))/log(2));
tau=2.^n;
l=2*pi;
x=0:h:2*pi;

T1=linspace(0,T,T/tau(1));
T2=linspace(0,T,T/tau(2));
T3=linspace(0,T,T/tau(3));

%% Classic fourier transform
%just a selector, so i can choose whether i'm solving for a=1, 10 or 100
apos=1;

switch apos
    
    case 1
    TT=T1;
    
    case 2
        TT=T2;
       
    case 3
        TT=T3;
end

du=@(x,xi,t) 1./(2*sqrt(pi*a(apos).*t)).*(sin(xi)+1/10.*sin(10*xi)).*exp(-(x-xi).^2./(4*a(apos).*t));

u=@(x,t)integral(@(xi) du(x,xi,t),-inf,inf);

aUvals=zeros(length(TT),length(x));
aUvals(1,:)=sin(x)+1/10*sin(10*x);
for i=2:length(TT)
    for j=1:length(x)

    aUvals(i,j)=u(x(j),TT(i));
    end
end

%% Fourier series

%solving only for a=1
  
n=10;
uk=zeros(length(T1),length(x));
U=uk;

for k=1:10
   I=(sin(1-k*pi/l)*l)/(2*(1-k*pi/l))-(sin(1+k*pi/l)*l)/(2*(1+k*pi/l));
   
   if (k==l/pi)
       I=1/2*l-1/4*sin(2*l);
   end
 
   II=1/10*((sin(10-k*pi/l)*l)/(2*(10-k*pi/l))-(sin(10+k*pi/l)*l)/(2*(10+k*pi/l)));
   
   if (k==10*l/pi)
      II=1/10*(1/2*l-1/40*sin(20*l));
   end
   
   III=2/l*(I+II);
   for i=1:length(T1)
       for j=1:length(x)
           
        uk(i,j)=III*sin(k*pi*x(j)/l)*exp(-(k*pi/l)^2*T1(i));
       end
   end
    U=U+uk;
end

%% FFT

xf=(-l/2):h:(l/2);
Nf=100;

y=(2*pi/l)*[-N/2:N/2];

y=fftshift(y);

fi=sin(xf)+1/10*sin(10*xf);

fiT=fft(fi);

uT=zeros(length(T1),length(x));
u=uT;
for i=1:length(t)
    uT(i,:)=exp(-a(1)*(y.^2)*T1(i)).*fiT;
    u(i,:)=ifft(uT(i,:));
end