clear all;

% The documentation of dst says:
% For best performance speed, the number of rows in x should be 2^m-1, for some integer m
N=2^8-1 % N=255

dx=pi/(N+1); % Grid spacing
x=linspace(0,pi,N+2);

u=zeros(N+2,1); % Solution of PDE

% Initial condition is a step function
u(1 : (N+1)/2)=0;
u((N+1)/2+1)=pi/2;
u((N+1)/2+2 : (N+2))=pi;

figure(1); clf;
plot(x,u,'k-');
xlim([0,pi]); ylim([0,pi]);
pause(5.0);

% Recall that dst does not take as inputs the values at x=0 and x=Pi since those are zero
c0=idst(u(2:N+1)); % Fourier coefficients of initial condition computed numerically

k2=([1:N].^2)'; % Eigenvalues of negative Laplacian with Dirichlet BCs

dt=0.001;
for timestep=1:20
   t=timestep*dt;

   % The vector of Fourier coefficients c is the solution in the basis
   % of eigenfunctions sin(n*x) which here become eigenvectors sin(pi*k*n/(N+1)) for k=1:N
   c = c0.*exp(-k2*t); % We can solve each ODE analytically here
   
   % Convert back to the standard basis in order to plot the solution:
   u(2:N+1) = dst(c); % Evaluate the solution on the grid of points (fast)
   u(1)=0; u(N+2)=0; % Impose the boundary conditions

   plot(x,u,'k-');
   xlim([0,pi]); ylim([0,pi]);
   pause(0.5);
end
