a=menu('Gitterpunkte in x-Richtung','10','20','40','80')-1;    
b=menu('Zeitschrittweite in s','1','2','4','8','16','32','64','128')-1;

% DER Befehl "menu3 erzeugt ein Fenster, wo man 10, 20, 40, 80 anclicken kann,
% diesen wird dann die Zahl 1,2,3 bzw. 4 zugeordnet


N=10*2^a;     % Anzahl der Gitterpunkte in Raumrichtung , d.h. x
M=4096/2^b;   % Anzahl der Gitterpunkte in Zeitrichtung, d.h. t 
h=1/N;        % Ortsschrittweite
tau=2^b;      % Zeitschrittweite
lambda=385;   
c=394;
rho=8930;     % Konstanten gemaess der DGL

% Setzen Anfangsbedingung und der Randbedingung
T=[];           
T=zeros(N+1,M+1);      % Ordnet T die Nullmatrix der Dimension N+1 x M+1 zu.
T(1:N+1,1)=100*sin(pi*[0:N]/N);

% Die letzte Zeile fuehrt eigentlich den Befehl
%
%  for i=1:N+1
%    T(i,1)=100*sin((i-1)*pi/N);
%  end
%
% aus, ist aber schneller als eine gewoehnliche for-Schleife


r=lambda/(c*rho)*tau/h^2;
for j=2:M+1   % Zeitschleife
    for i=2:N  % ueber alle inneren Ortspunkte
      T(i,j)=(1-2*r)*T(i,j-1)+r*(T(i-1,j-1)+T(i+1,j-1));  % Rekursionsformel gemaess Vorlesung
      
    end    
end    


% Visualisierung mittels plot-Befehl, siehe "help plot"
% Es wird T nur zu gewissen Zeitpunkten gezeichnet.
s=(M)/8;
x=[0:N]/N;
plot(x,T(1:N+1,0*s+1),'k',x,T(1:N+1,s+1),'b',x,T(1:N+1,2*s+1),'g',...
    x,T(1:N+1,3*s+1),'r',x,T(1:N+1,4*s+1),'c',x,T(1:N+1,5*s+1),'y',...
    x,T(1:N+1,6*s+1),'m',x,T(1:N+1,7*s+1),'k',x,T(1:N+1,8*s+1),'b')
legend('0 s','512 s','1024 s','1536 s','2048 s','2560 s','3072 s','3584 s','4096 s')
