% Waermeleitungsgleichung im Zylinder - rotationsymmetrisch
% 
% - explizite Methode -
%
% r=1 : Dirichlet (=fr1!)
% r=0 : Neumann (=0!) - Symmetrie bzgl. r=0!!
% t=0 : Anfangsbedingung (=fr!)
%
% h : Maschenweite fuer r
% k : Maschenweite fuer t
% n0 : Anzahl der Iterationen
%
% R. MOHR 3. 3. 1998
%
h=0.025;k=0.0001; % Maschenweite
n0=6000; % Iterationen
n00=3850; % Beobachtungszeitpunkt
fr='zeros(size(x))'; % Anfangsbedingung
fr1='exp(-t)-1'; % Dirichlet r=1
%fr1='sin(4*t)';
x=[0:h:1];
n=length(x);
eval(['xi=' fr ';']);
t0=n0*k;
[xg,yg]=meshgrid([0:h:1],[0:100*k:t0]);
l=0;
f1=1-2*k/(h*h);
f2=k/(h*h);
f3=k*(n-1)/(2*h);
zg=[];
xii=zeros(size(xi));
for iz=1:(n0+1),
  t=(iz-1)*k;
  xii(1)=f1*xi(1)+(1-f1)*xi(2);
  eval(['xii(n)=' fr1 ';']);
  for izz=2:(n-1),
    xii(izz)=f1*xi(izz)+(f2+f3/(izz-1))*xi(izz+1)+(f2-f3/(izz-1))*xi(izz-1);
  end;
  xi=xii;
  if rem(l,100)==0, zg=[zg;xi];end;
  l=l+1;
end;

mesh(xg,yg,-zg);title(['Waermeleitungsgleichung im Zylinder  fr= ' fr1 ';  g= ' fr ';    k= ' num2str(k) '; h= ' num2str(h)]); 
xlabel('r');ylabel('t');zlabel('U(r,t)');
for jz=[8:8],
clear u uu rr;
nt=round(n00/100)+jz*2;
t0=nt*100*k;

uu=zg(nt,1:n);rr=[0:h:1];
[x,y]=meshgrid([-1.25:h:1.25]);
[n1,n2]=size(x);
for iz=[1:n1],
  for izz=[1:n2],
    r=sqrt(x(iz,izz).^2+y(iz,izz).^2);
    if r>1, u(iz,izz)=uu(n);,
    else,
      ir=find((rr>=r)&(rr<r+h));
      ir=min(ir);
      %u(iz,izz)=uu(ir);
    end;
  end;

end;
%pause;close;
figure;
mesh(x,y,u);
xlabel('x');ylabel('y');zlabel('U(x,y,t0)');
title(['Waermeleitung im Zylinder  f= ' fr1 ';  t0= ' num2str(t0) ';    k= ' num2str(k) '; h= ' num2str(h)]); 
end;