L = pi;
T = 20;
n_t = 10000;
n_x = 100;
x = linspace(0, L, n_x);
t = linspace(0, T, n_t);

delta_t = t(2) - t(1);
delta_x = x(2) - x(1);

A_forward = zeros(n_x, n_x );
A_backward = zeros(n_x, n_x);

lambda = delta_t / (delta_x)^2;

for i = 1 : n_x 
    A_forward(i,i) = 1 - 2 * lambda;
    A_backward(i,i) = 1 + 2 * lambda;
    if (i-1 > 0) && (i+1 < n_x+1)    
        A_forward(i,i-1) = lambda;
        A_forward(i,i+1) = lambda;
        A_backward(i,i-1) = -lambda;
        A_backward(i,i+1) = -lambda;
    end
end

% Neumann condition
A_forward(1,2) = 2 * lambda;
A_backward(1,1) = 1- lambda;
A_forward(n_x,n_x-1) = 2 *  lambda;
A_backward(n_x,n_x-1) = 1-  lambda;

A_inv_backward = inv(A_backward);
u_forward = zeros( n_x, n_t);
u_backward = zeros( n_x, n_t);

%Initialize u(x) = sin(x)
u_forward(:,1) = sin(x);
u_backward(:,1) = sin(x);

for i = 2 : n_t
    u_forward(:,i) = A_forward * u_forward(:,i-1);
    u_backward(:,i) = A_inv_backward * u_backward(:,i-1);
end

plot(u_backward)
lambda
% % 
% % plot(u_backward(:,5));
% % hold on;
% % plot(u_backward(:,20));
% % plot(u_backward(:,50));
% % plot(u_backward(:,100));
% % plot(u_backward(:,500));
% % plot(u_backward(:,1000));
