close all; clear all;

L=1; k=0.001; c=1; rho=1; tfinal =150;


%n=10; 
n=input('Enter n=')
h=(L-0)/n;
for i=1:n+1,
   x(i)=0 + (i-1)*h;
end
for i=1:n+1                  % Initial Condition
   u0(i) = sin(pi*x(i));
end
a=rho*c; b=k/a;

plot(x,u0)
pause(3)
hold

dt = a*h*h/(2.1*k); kmax=fix(tfinal/dt);

for kt=1:kmax,              % Time marching
   u1(1) = 0; u1(n+1)=0;     % BC
   for i=2:n
      u1(i) = dt/a + u0(i)*(1-2*dt*b/(h*h)) + ...
         (dt*b/(h*h))*(u0(i-1) + u0(i+1));
   end
   u0=u1;
   plot(x,u1)
   pause(2)
end
plot(x,u1)