
clear

X = linspace(0,1,101);
T = linspace(0,3,91);

u0 = @(x) max(1-8*abs(x-.3),0);
v0 = @(x) 0;

m = 20;
for n=1:m
    a(n) = 2*integral(@(x) u0(x).*sin(n*pi*x),0,1);
    b(n) = 2/(n*pi)*integral(@(x) v0(x).*sin(n*pi*x),0,1);
end
for i=1:length(T)
    U(:,i) = sin(n*pi*X).*(a*cos((1:m)'*pi*T(i))+b*sin((1:m)'*pi*T(i)));
end

figure(1)
plot(1:length(a),a,'-x','linewidth',1,'markersize',12);
axis tight
xlabel('$x$','interpreter','latex')
title('First 20 Fourier coefficients','FontWeight', 'Normal')
set(gca,'fontsize',18)

figure(2)
for i=1:20
    v(:,i) = a(i)*sin(i*pi*X);
end
plot(X,u0(X),'k',X,v*ones(m,1),X,v(:,1:5)*ones(5,1),'linewidth',1);
axis tight
xlabel('$x$','interpreter','latex')
title('Initial data and approximation','FontWeight', 'Normal')
l = legend('Initial data','Appr. $n=20$','Appr. $n=3$');
set(l,'interpreter','latex','location','northeast')
set(gca,'fontsize',18)

figure(3)
id = [1 6 16];
for i=id
    plot(X,v*cos((1:m)'*pi*T(i)),'linewidth',1);
    hold on
end
hold off
l = legend('$t=0$',['$t=$',num2str(T(id(2)),3)],['$t=$',num2str(T(id(3)))]);
set(l,'interpreter','latex','location','northeast')
xlabel('$x$','interpreter','latex')
ylabel('$u$','interpreter','latex')
title('Wave expansion & reflection','FontWeight', 'Normal')
set(gca,'fontsize',18)

figure(4)
for i=1:length(T)
    Z(:,i) = v*cos((1:m)'*pi*T(i));
end
mesh(X,T,Z')
axis tight
xlabel('$x$','interpreter','latex')
ylabel('$t$','interpreter','latex')
title('Solution wave equation','FontWeight', 'Normal')
set(gca,'fontsize',18)

figure(5)
for i=1:length(T)
    plot(X,Z(:,i),'k','linewidth',1)
    axis tight
    set(gca,'fontsize',18,'xlim',[0 1],'ylim',[-1 1])
    xlabel('$x$','interpreter','latex')
    ylabel('$u$','interpreter','latex')
    M(i) = getframe;
end

movie(M)
