clc; clear; clf
%% Input parameter
Nx = 20;
x = linspace(0,1,Nx);
h = x(2) - x(1);
alpha = 0.45;
k = alpha*h^2;
Nt = 100;
%% Initialization
u(1:Nx,1:Nt+1)=0;
u(:,1) = sin(pi*x);
%% Computation
for n=1:Nt
for i=2:Nx-1
u(i,n+1) = alpha*(u(i+1,n) - 2*u(i,n) + u(i-1,n)) + u(i,n);
end
end
figure(1)
plot(x,u(:,1),'ko-',x,u(:,100),'ro-')