function [] = dftcs4()
%DFTCS4 Diffusionsgleichung in 2D mit FTCS
%       Neumann Randbedingungen
%       u_x = 0  bzw. u_y = 0 am Rand fuer alle t

clc
close all
help dftcs4

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 = -L/2 bis x = +L/2
dx = L/(N-1);             % Schrittweite im Ortsgitter
dy = L/(N-1);             % Schrittweite im 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
y = (0:N-1)*dy - L/2;     % x(1) ,..., x(N)  : raeumliche Gitterpunkte

u    = zeros(N, N);       % u(1,1) ,..., u(1,N)  : Naeherungsloesung
                          % u(N,1) ,..., u(N,N)
uneu = zeros(N, N);


nstep = 1600;

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

% Plot
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)])

%% Konventionell (Doppelschleife)
%for istep = 1:nstep
%    for j = 2:(N-1)
%        for k = 2:(N-1)
%            uneu(j,k) = u(j,k) + ...
%                        rx * ( u(j+1,k  ) - 2*u(j,k) + u(j-1,k  ) ) + ...
%                        ry * ( u(j  ,k+1) - 2*u(j,k) + u(j  ,k-1) ) ;
%        end
%    end
%    u = uneu;
%    t = t + dt;
%end

% Das Gleiche mit Matrizen
for istep = 1:nstep
    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)) ) ;
    uneu(:,1) = u(:,1) + 2.0 * rx * ( u(:,2)   - u(:,1) );
    uneu(:,N) = u(:,N) + 2.0 * rx * ( u(:,N-1) - u(:,N) );
    uneu(1,:) = u(1,:) + 2.0 * ry * ( u(2,:)   - u(1,:) );
    uneu(N,:) = u(N,:) + 2.0 * ry * ( u(N-1,:) - u(N,:) );
    u = uneu;
    t = t + dt;
    if ( istep == nstep/4 )
       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 2])
xlabel('x')
ylabel('y')
zlabel('u(x,y)')
title(['2D-Diffusionsgleichung mit FTCS, t = ', num2str(t)])
colorbar

% Postscript
epsfile = 'dftcs4.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
%----------------------------------------------------------------------