function dufort;

clf;


tt=0.04;
h=0.01; plot([0:h:1],ftcs(h,tt),'b.'); hold on;
%h=0.02; plot([0:h:1],ftcs(h,tt),'r.'); hold on;
%h=0.05; plot([0:h:1],ftcs(h,tt),'g.'); hold on;
h=0.01; plot([0:h:1],df(h,tt),'r'); hold on;



function y=ftcs(h,tt);

n=round(1/h);
x=zeros(1,n);
x=sin([0:n]*2*pi/n);
%x(round(0.2/h))=10;
y=x;

% plot([h:h:1],y,'g'); hold on;

m=0.4;
k=m*h^2;
t=0;
while k*t<tt;
  for i=2:n
     y(i)=x(i)+m*(x(i-1)-2*x(i)+x(i+1)); 
  end
  
%  plot([0:h:1],y,'b'); axis([0 1 -1 1]); pause;
  
  x=y;
  t=t+1;
end


function z=df(h,tt);

n=round(1/h);
x=zeros(1,n);
x=sin([0:n]*2*pi/n);
y=x;
z=x;
m=0.4;
k=m*h^2;
k=m*h^2*20;

m=k/h^2; 

for i=2:n
   y(i)=x(i)+m*(x(i-1)-2*x(i)+x(i+1)); 
end

t=0;

% plot([h:h:1],y,'g'); hold on;

%hold off;
while k*(t+1)<tt;
%  plot([0:h:1],z,'w','linewidth',2);

  for i=2:n
     z(i)=((1-2*m)*x(i)+2*m*(y(i-1)+y(i+1)))/(1+2*m); 
  end
  
%  plot([0:h:1],z,'r','linewidth',2); axis([0 1 -1 1]); pause;
  x=y;
  y=z;
  t=t+1;
end