%% Method of lines for the heat equation


%% Parameters

d = 1;                  % the length of the spatial interval
n = 20;                 % the number of points on the spatial interval
h = d/n;                % the length of a spatial segment

a = 1;

%% Building the state-space matrices

c = zeros(n-1,1);
r = zeros(n-1,1);

c(1:2) = [-2 1];
r(1:2) = [-2 1];

A = toeplitz(c,r);
A = a/h^2 * A;

B = zeros(n-1,2);
B(1,1) = 1;
B(n-1,2) = 1;
B = a/h^2 * B;

C = eye(n-1,n-1);
D = zeros(n-1,2);

G = ss(A,B,C,D);

%% Analyzing the spectrum of the state space matrix 

figure(1)
eA = eig(A);
plot(real(eA),imag(eA),'xr')
xlabel('Real')
ylabel('Imaginary')

%% Examining the sparsity structure of the matrix A

figure(2)
spy(A)
set(gca,'FontSize',18)

% print -depsc2 ../../lectures/figures/spy_A.eps
% !epstopdf ../../lectures/figures/spy_A.eps

%% Simulation

T0 = zeros(n-1,1);
% T0([4 5 6],1) = 10;                     % the initial temperature profile
% TO((n/2-n/10):(n/2+n/10),1) = 100

tf = 0.5;
t=linspace(0,tf,100);

Ta = 10;                                % the temperature at the left end
Tb = 0;                                 % the temperature at the right end
Tab = repmat([Ta;Tb],1,length(t));      

[T,t,x] = lsim(G,Tab,t,T0);

F(length(t)) = struct('cdata',[],'colormap',[]);

figure(3)
for it = 1:length(t)
    plot(0:h:d,[Ta, T(it,:), Tb],'.-',[0,d],[Ta,Tb],'r.','MarkerSize',10)
    xlabel('x')
    ylabel('\theta(x,t)')
    axis([-1,d+1,0,10])
    grid on
    text(-0.5,7,[' t = ',num2str(t(it))],...
            'VerticalAlignment','middle',...
            'HorizontalAlignment','left',...
            'FontSize',12)
    F(it) = getframe(gcf);
end

%% Snapshots at a few time instances
figure(4)
subplot(1,3,1)
    it = 1;
    plot(0:h:d,[Ta, T(it,:), Tb],'.-',[0,d],[Ta,Tb],'r.','MarkerSize',10)
    set(gca,'FontSize',14)
    xlabel('x')
    ylabel('\theta(x,t)')
    axis([0,d,0,10])
    grid on
    text(0.5,7,[' t = ',num2str(t(it),2)],...
            'VerticalAlignment','middle',...
            'HorizontalAlignment','left',...
            'FontSize',14)
subplot(1,3,2)
    it = 10;
    plot(0:h:d,[Ta, T(it,:), Tb],'.-',[0,d],[Ta,Tb],'r.','MarkerSize',10)
    set(gca,'FontSize',14)
    xlabel('x')
    ylabel('\theta(x,t)')
    axis([0,d,0,10])
    grid on
    text(0.5,7,[' t = ',num2str(t(it),2)],...
            'VerticalAlignment','middle',...
            'HorizontalAlignment','left',...
            'FontSize',14)
subplot(1,3,3)
    it = 100;
    plot(0:h:d,[Ta, T(it,:), Tb],'.-',[0,d],[Ta,Tb],'r.','MarkerSize',10)
    set(gca,'FontSize',14)
    xlabel('x')
    ylabel('\theta(x,t)')
    axis([0,d,0,10])
    grid on
    text(0.5,7,[' t = ',num2str(t(it),2)],...
            'VerticalAlignment','middle',...
            'HorizontalAlignment','left',...
            'FontSize',14)        

% print -depsc2 ../../lectures/figures/mol_3_responses_BC.eps
% !epstopdf ../../lectures/figures/mol_3_responses_BC.eps

%% Playing the animation
figure(5)
movie(gcf,F,3,10)