function solveheate=solveheate(u,t);
%solves the heat equation u_t=u_xx on (0,\pi) with zero boundary condition.
% (for solving u_t=k u_{xx} on (0,L) one just needs to rescale time; it is a good exercise to work out how exactly)
%
%input variables:
%
%    u   field of the initial datum (temperature at time t=0 at points 1,2,...,N)
%    t   time at which we wish to calculate the temperature
%
% output variables:
%
% temperature field at time t (same size as the input)
%
n=size(u,2); %for optimal performance of fft (which may not be our main concern in 1d problems
              %with the heat equation) it is best to take n=2^m;
u(n+1:2*n)=-u(1:n); %extension to an odd function to accommodate the zero boundary condition
N=[0:n-1]; %start setting up the multiplier of the Fourier modes which will be used for the solution
           %this part takes care of what would be ``positive frequencies"
           %in full Fourier series; in the discrete case the frequancies
           %are defined only modulo 2n
N(n+1:2*n)=-flip([1:n]); %this part deals with the ``negative frequencies"
                         %the particular ordering has to do with how exactly the fft is 
                         %defined in matlab
                         %If we were working with a more sensitive equation
                         %at the limit of the grid resolution, we might
                         %consider setting N(n)=0. If we leave N as it is, it is
                         %slighly non-symmetric: N(n-1)=n-1, N(n)=-n,
                         %N(n+1)=-(n-1); setting N(n)=0 looks better for
                         %many reasons. However, for the heat equation it
                         %does not matter
M=exp(-t*(N.*N)); %the multiplier which captures the time evolution on the Fouriet side
fu=fft(u); %take the DFT; 
fut=M.*fu; %solution at time t on the Fourier side
ut=ifft(fut); %invert the DFT
solveheate=ut(1:n); %the relevant part which gives the solution at time t
end


