close all; clear all
L=1; vel=0.1; dec=0.001; tfinal=15;

n = input('n= ')

h = L/n;
for i=1:n+1
   x(i) = (i-1)*h;
   %u0(i) = (i<(n/2+1))*sin(pi*x(i)*2)+(i>n/2)*0;
   u0(i)=cos(pi*x(i));
end
plot(x,u0); hold
dt = h/(vel*1.1+h*dec*1.1);
kmax=fix(tfinal/dt); a=dt/h; t=0;
for k=1:kmax
   u1(1)=exp(-t);
   for i=2:n+1
      u1(i)=a*vel*u0(i-1)+(1-a*vel-dec*dt)*u0(i)...
         +dt*my_f(x(i),t,vel,dec);
   end
   u0=u1; t=t+dt;
   plot(x,u1); pause(0.5)
   %waitforbuttonpress
end
for i=1:n+1
   err(i)=u1(i)-u1exact(x(i),t);
end
figure(2); plot(x,err)
e2=[norm(err,1),norm(err,2),norm(err,inf)]