% Waermeleitungsgleichung mit Dirichlet-RBen
% - implizite Euler-Methode -
%
% x=0 : Dirichlet        (=fr1!)
% t=0 : Anfangsbedingung (=fr!)
%
% h : Ortsschrittweite
% k : Zeitschrittweite
% n0 : Anzahl der Iterationen*10
%
%
h=0.05;k=0.02; % Gitterweiten
n0=floor(1/k); % Iterationen

fr = @(x) ones(size(x));                % Anfangsbedingung
f = @(t,x) (-1)*(t.^2+x).*(exp(-t.*x)); % inhomogene rechte Seite
fr1 = @(t) exp(-t);                     

x=[0:0.05:1];n=length(x);
xi = fr(x);
xi = xi';

n0=10*n0;t0=n0*k;
[xg,yg]=meshgrid([0:h:1],[0:k:t0]);

%-- Aufstellung der Steifigkeitsmatrix
r=k/(h^2);

a=zeros(n,1);
b=zeros(n,1);
c=zeros(n,1);

a(1)=1;
a(2:n)=2+2*r;

b(3:n)=-r;

c(1:n-2)=-r;
c(n-1)=-2*r;


A=spdiags([c a b], -1:1, n, n);


%-- Systemmatrix fŸr impliziten Euler
B=-A;
for iz=1:n,
  B(iz,iz)=2-2*r;
end;


cpe=cputime;
l=0;
zg = zeros(n0+1,size(xi,1));
for iz=1:(n0+1),
    t=(iz-1)*k;

    fkk=f(t,x);
    xr1 = fr1(t);
    t=iz*k;
    fkk1 = f(t,x);

    xr11 = fr1(t);
    d=B*xi;
    d=d-k*(fkk'+fkk1');

    d(1)=xr11;
    xi=A\d;
    zg(iz,:)=xi';
end;
ept=cputime-cpe


mesh(xg,yg,zg); 
xlabel('x');ylabel('t');zlabel('U(x,t)');
