clear;
h = 0.1;
T = 32*pi;
N = ceil(T/h);
M = 15;

A11 = zeros(M,M);
A12 = eye(M);
A21 = zeros(M,M);
for j = 1:M
    A21(j,j) = -2;
end
for j = 2:M
    A21(j-1,j) = 1;
    A21(j,j-1) = 1;
end
A22 = zeros(M,M);

A = [A11 A12; A21 A22];

z = [ (M+1)/2-abs( (1:M)-(M+1)/2 )' ; zeros(M,1)];

lim = [0 M+1 -max(z) max(z)];

plot(0:M+1,[0; z(1:M); 0]);
axis(lim);
getframe;

for j = 1:N
    
    %z = z + h*A*z;
    z = (eye(2*M)-h*A)\z;
    
    plot(0:M+1,[0; z(1:M); 0]);
    axis(lim);
    getframe;
    
end

