% This is an example to understand the methods proposed in the book in section 6.1
% programmer: C. Vuik
% date: 26-03-2020
% e-mail: c.vuik@tudelft.nl

clear
clf

n = 100;
h = 1/n;

% possible choices for the time step

% dt = 1.1/n;
dt = 1/n;
% dt = 0.99/n;
dt = 0.8/n;
teind = 5;

for j = 1:n
x(j) = (j-1)*h;
a(j+1,j) = 1/(2*h)+dt/(2*h*h);
a(j,j) = -dt/(h*h);
a(j,j+1) = -1/(2*h)+dt/(2*h*h);
end

x(n+1) = j*h;
a(n+1,n+1) = -dt/(h*h);
a(1,n+1) = 1/(2*h)+dt/(2*h*h);
a(n+1,1) = -1/(2*h)+dt/(2*h*h);

q(n+1,1) = 0;

shift = floor(n/8);

for j = n/4:n/2

q(j-shift,1) = (1+sin(8*pi*j*h-pi/2))/2;
q(j+n/2-shift,1) = 1;
end

plot(x,q);
hold on

% the time loop

for j = 1:floor(teind/dt)

q = q+dt*a*q;

end

plot(x,q,'o');