function [] = dftcs1()
%DFTCS1  Diffusionsgleichung in 1D mit FTCS
%        AB: Rechtecksverteilung der Breite delta
%        RB: Dirichlet Randbedingungen
%            u(-L/2,t) = u(+L/2,t) = 0

clc
help dftcs1

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
D  = 1.0;                 % Diffusionskoeffizient

r = D*dt/dx^2;

if ( r > 0.5 )
   disp(' WARNUNG: FTCS ist instabil')
end

x = (0:N-1)*dx - L/2;     % x(1) ,..., x(N)  : raeumliche Gitterpunkte
                          % u(1) ,..., u(N)  : Naeherungsloesung

nstep = 300;
nplots = 30;
plot_step = nstep/nplots; % Anzahl der Zeitschritte zwischen den Plots

% Anfangsbedingung
u = f(x);                 % u(x,0)= f(x)
iplot = 1;
tplot(iplot) = 0.0;
uplot(iplot,:) = u;
iplot = iplot + 1;
uneu = zeros(1, N);

for istep = 1:nstep
    % Konventionell
    % for j = 2:(N-1)
    %     uneu(j) = u(j) + r * (u(j+1) - 2.0*u(j) + u(j-1));
    % end
    uneu(2:(N-1)) = u(2:(N-1))  + ...
                    r * (u(3:N) - 2*u(2:(N-1)) + u(1:(N-2)));
    u = uneu;
    if ( mod(istep,plot_step) == 0 )
       uplot(iplot,:) = u;
       tplot(iplot) = istep*dt;
       iplot = iplot + 1;
    end
end

% Plot
figure(1)
mesh( x, tplot, uplot )
view(40,20)
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
title('1D-Diffusionsgleichung mit FTCS')

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

%----------------------------------------------------------------------
function [u] = f(x)
% Anfangsbedingung
delta = 0.1;
N = length(x);
u = zeros(1,N);
for i = 1:length(x)
    if ((x(i) < -delta/2) || (x(i) > delta/2))
         u(i) = 0.0;
    else
         u(i) = 1.0/delta;
    end
end
%----------------------------------------------------------------------
