function [] = dftcs3()
%DFTCS3  Diffusionsgleichung in 2D mit FTCS
%        Dirichlet Randbedingungen
%        u(x,y,t) = 0   am Rand fuer alle t

clc
help dftcs3

N  = input(' Anzahl der raeumlichen Gitterpunkte N [40]: ');
dt = input(' Zeitschritt dt [1.0e-4]: ');

if isempty(N)
   N = 40;
end
if isempty(dt)
   dt = 1.0e-4;
end
L  = 1.0;                 % x,y = -L/2 bis x,y = +L/2
dx = L/(N-1);             % Schrittweite im x-Ortsgitter
dy = L/(N-1);             % Schrittweite im y-Ortsgitter
D  = 1.0;                 % Diffusionskoeffizient
rx = D*dt/dx^2;
ry = D*dt/dy^2;

if ( (rx + ry) > 0.5 )
   disp(' WARNUNG: FTCS ist instabil')
end

x = (0:N-1)*dx - L/2;     % x(1) ,..., x(N)     : raeumliche Gitterpunkte x
y = (0:N-1)*dy - L/2;     % y(1) ,..., y(N)     : raeumliche Gitterpunkte y
u    = zeros(N, N);       % u(1,1) ,..., u(1,N) : Naeherungsloesung
                          % u(N,1) ,..., u(N,N)
uneu = zeros(N, N);
nstep= 200;               % nstep Zeitschritte

% Anfangsbedingung
t = 0.0;
u = f(x,y,dx);            % u(x,y,0)= f(x,y)

% Graphik
figure(1), mesh( x, y, u ), view(40,30)
xlabel('x'), ylabel('y'), zlabel('u(x,y)')
title(['2D-Diffusionsgleichung mit FTCS, t = ', num2str(t)])

for istep = 1:nstep
    % Konventionell (Doppelschleife)
    % for i = 2:(N-1)
    %     for j = 2:(N-1)
    %         uneu(i,j) = u(i,j) + ...
    %                     rx * ( u(i+1,j  ) - 2*u(i,j) + u(i-1,j  ) ) + ...
    %                     ry * ( u(i  ,j+1) - 2*u(i,j) + u(i  ,j-1) ) ;
    %     end
    % end
    uneu(2:(N-1),2:(N-1)) = u(2:(N-1),2:(N-1)) + ...
    rx*( u(2:(N-1),3:N) - 2*u(2:(N-1),2:(N-1)) + u(2:(N-1),1:(N-2)) ) + ...
    ry*( u(3:N,2:(N-1)) - 2*u(2:(N-1),2:(N-1)) + u(1:(N-2),2:(N-1)) ) ;
    u = uneu;
    t = t + dt;
    if ( istep == nstep/2 )
       figure(2)
       mesh( x, y, u )
       view(40,30)
       xlabel('x'), ylabel('y'), zlabel('u(x,y)')
       title(['2D-Diffusionsgleichung mit FTCS, t = ', num2str(t)])
    end
end

% Graphik
figure(3), mesh( x, y, u )
view(40,30), axis([-L/2 L/2 -L/2 L/2 0 8])
xlabel('x'), ylabel('y'), zlabel('u(x,y)')
title(['2D-Diffusionsgleichung mit FTCS, t = ', num2str(t)])

% Postscript
epsfile = 'dftcs3.eps';
print('-depsc', epsfile)

%----------------------------------------------------------------------
function [u] = f(x,y,dx)
% Anfangsbedingung
delta = 2.0*dx; u = zeros(length(x),length(y));
for i = 1:length(x)
    for j = 1:length(y)
        if ( ((x(i) >= -delta/2) && (x(i) <= delta/2)) && ...
             ((y(j) >= -delta/2) && (y(j) <= delta/2)) )
             u(i,j) = 1.0/delta^2;
        else
             u(i,j) = 0.0;
        end
    end
end
%----------------------------------------------------------------------