% This code solve the heat equation on a long thin rod (L=length )
% u_t = a^2 u_xx
% u(x,0) = f(x)
% u(0,t) = u(L,t)
clear all

a = 1;
L = 25;
f = @(x,n,L) -x.*(x-25).*sin(n*pi*x/L);

% g = @(x)  ;
% 10*sin(pi*x).;
% number of terms in the partial sum
N = 50;
% how far time
tmax = 100;
% time step
delt = 1;
% space step
delx = 0.5;

x = [0:delx:L];
t =[0:delt:tmax];
Nt = length(t);

% compute fourier coefficients
A = zeros(N,1);
for n = 1:N
    A(n) = 2/L*quad(f,0,L,[],[],n,L);
end
% q = quad(fun,a,b,tol,trace)

%compute heat
u = zeros(Nt,length(x));
%sweep over time
for j = 1:Nt
% sweep over space
    for k= 1:length(x)
% computation of heat at a given place x(k) at time t(j)
        for n = 1:N
            mu = n*pi/L;
            u(j,k) = u(j,k) + A(n)*exp(-a^2*mu^2*t(j))*sin(mu*x(k));
%             u(j,k) = u(j,k) + A(n)*sin(mu*x(k));

        end
    end
end


for j = 1:Nt
    plot(x,u(j,:))
    pause
end


% [X,T] = meshgrid(x,t);
% mesh(X,T,u)