clc;
close all;
clear all;

n = 19;
L = 400;
beta = 1;

dx = (2*pi)/(n+1);
xmaxim = 2*pi;
m = (xmaxim/dx)-1;

dt = 10/L;
tmaxim = 10;

lamda = (beta*dt)/(dx^2);
c0 = 1+2*lamda;
T = zeros(m,1);
j = 1;

for x = dx:dx:(xmaxim-dx)
    T(j,:) = sin(x);
    j = j+1;
end

A = zeros(m,m);
for i = 2:m-1
    A(i,i)   = c0;
    A(i,i-1) = -1*lamda;
    A(i,i+1) = -1*lamda;
end
A(1,1)  = c0;
A(m,m)  = c0;
A(1,2)  = -1*lamda;
A(m,m-1)= -1*lamda;
B = zeros(tmaxim/dt+1, m+2);
k=2;
B(1,:) = [0 T.' 0];
for t = dt:dt:tmaxim
    T1 = A\T;
    B(k,:) = [0 T.' 0];
    T = T1;
    k = k+1;
end
X = 0:dx:xmaxim;
Y = 0:dt:tmaxim;
%mesh (X:Y:B);
figure
subplot(3,1,1)
plot(Y,B(:,2));
AXIS([0 10 0 1]);

subplot(3,1,2)
plot(Y,B(:,5));
AXIS([0 10 0 1]);

subplot(3,1,3)
plot(Y,B(:,10));
AXIS([0 10 0 1]);


%Tdeltax1    = B(:,2);
%Tdeltax4    = B(:,5);
%Tdeltax10   = B(:,11);