% solves string problem,
% u_tt - c^2*u_xx = 0 for 0<x<1, t>0,
% u(0,t) = u(1,t) = 0 for t >= 0,
% u(x,0) = u0(x) for 0 <= x <= 1,
% using truncated Fourier series

c=1; 
N = 100; % number of terms in truncated Fourier series
numx = 8*N; % number of x values
numt = 300; % number of t values
nvec = (1:N); 
tvec = linspace(0,6/c,numt)'; % time lenght = three periods
x = linspace(0,1,numx)';

%%% initial position of string
u0 = (x.*(1-x).^2); 

%%% solve with truncated Fourier series
Phi = sin(pi*x*nvec);
T = cos(pi*c*tvec*nvec);
UU = (Phi*diag(Phi'*u0)*T')'*2/(numx-1);

%%% plot solution u(x,t) as a surface
figure(1), clf
%mesh(x,tvec,UU)
mesh(linspace(0,1,N),tvec(1:2:end),UU(1:2:end,1:8:end))
xlabel('x','fontsize',16)
ylabel('t','fontsize',16)
zlabel('u(x,t)','fontsize',16)
title(sprintf('string solution when c = %f',c),'fontsize',16)

%figure(3), clf
%plt = waterfall(linspace(0,1,N),tvec(1:4:end),UU(1:4:end,1:8:end));
%set(plt,'edgecolor','b')     
%xlabel('x','fontsize',16)
%ylabel('t','fontsize',16)
%zlabel('u(x,t)','fontsize',16)
%ylim([min(tvec) max(tvec)])
%title(sprintf('string solution when c = %f',c),'fontsize',16)


%%% animated plot of string position
figure(2), clf
ht = plot(x,UU(1,:));
axis([0 1 -1.2*max(abs(u0)) 1.2*max(abs(u0))]) 
set(gca,'DataAspectRatio',[1 1 1])
xlabel('x','fontsize',16)
ylabel('t','fontsize',16)
title(sprintf('t = %5.3f',tvec(1)),'fontsize',16)
input('hit return to start\n');
for j = 2:numt
    set(ht,'YData',UU(j,:))
    title(sprintf('t = %5.3f',tvec(j)),'fontsize',16)
    drawnow
%    pause(0.03)
end