%% Definition des Gitters 
clc; clear; close all;
                        %Grenzen des Gitters 
LR=0;                   %Begrenzung links 
RR=1;                   %Begrenzung rechts

M=40;                   %Vorgabe der Schritte in x-Richtung
 
x=linspace(LR,RR,M+1);  %Ermittlung der Schrittweite in x-Richtung
dx=x(2)-x(1);

K=500;                  %Vorgabe der Endbedingung K

dt=0.025*dx^2;            %Ermittlung der Schrittweite in t-Richtung
                        %bei konstantem dx ist Fehler umso größer, desto größer dt ist
                        %Stabilitätsbedingung explizit laut Numerische Mathematik von
                        %Schwarz S. 453 Gleichung (10.50)
                        %dt=0.5*dx^2;

t_end=dt*K;             %Ermittlung des t-Vektors
t=0:dt:t_end;



aa=zeros(K+1,M+1);      %Definition der Matrix a(x,t)
r=zeros(K+1,M+1);       %Definition der Matrix r(x,t)
f=zeros(K+1,M+1);       %Definition der Matrix f(x,t)

%% Variantenwahl
Vari=1;                                     %|| oder Abfrage

if Vari==1                                  %homogen, a=konstant, r=konstant
    for jj=1:K+1
        for ii=1:M+1
            aa(jj,ii)=1;                    %konstante Temperaturleitzahl
            r(jj,ii)=aa(jj,ii)*dt/dx^2;
            f(jj,ii)=0;
        end
    end
elseif Vari==2                              %inhomogen, a=konstant, r=konstant
    for jj=1:K+1
        for ii=1:M+1
            aa(jj,ii)=1;                    %konstante Temperaturleitzahl
            r(jj,ii)=aa(jj,ii)*dt/dx^2;
            f(jj,ii)=500*x(ii)+1*t(jj);
        end
    end
elseif Vari==3                              %inhomogen, a=variabel, r=variabel
    for jj=1:K+1
        for ii=1:M+1
            aa(jj,ii)=0.05*t(jj)+0.05*x(ii);%veränderliche Temperaturleitzahl abhängig von t und x
            r(jj,ii)=aa(jj,ii)*dt/dx^2;
            f(jj,ii)=500*x(ii)+1*t(jj);
        end
    end
elseif Vari==4                              %Offen für ein freies Beispiel
    u_ana=@(x,y) x.^2+y.^2;
end

%% FTCS

%Start- und Randbedingungen festlegen
    uIC=sin(pi*x);         %Anfangsbedingung
    cc=1*t;                %linker Rand: Neumannrand-Bedingung                 
    uR=0.2*t;              %rechter Rand: Dirichlet

if r<=0.5
    
    U1=zeros(K+1,M+1);
    U1(1,:)=uIC;            %IC
    U1(:,M+1)=uR;           %BC

    for jj=2:K+1
        for ii=M:-1:2
        % FTCS inhomogen: U1(jj,ii)=r*U_W + (1-2*r)*U_P + r*U_E + dt*f(jj,ii)
        U1(jj,ii)=r(jj-1,ii-1)*U1(jj-1,ii-1)+(1-2*r(jj-1,ii))*U1(jj-1,ii)+r(jj-1,ii+1)*U1(jj-1,ii+1)+dt*f(jj,ii);
        end
        % Neumann-Bedingung inhomogen: U1(jj,1)=2*r*U_E+(1-2*r)*U_P+r*2*dx*cc(t)+dt*f(jj,ii) 
        U1(jj,1)=2*r(jj-1,2)*U1(jj-1,2)+r(jj-1,2)*2*dx*cc(jj-1)+(1-2*r(jj-1,1))*U1(jj-1,1)+dt*f(jj,ii);
    end
else
    disp('Für das FTCS-Verfahren sind die Stabilitätskriterien nicht erfüllt');
end

%% BTCS Verfahren 

a=zeros(K,M);       %Matrix für a
for jj=1:K
    for ii=2:M
        a(jj,ii)=-r(jj+1,ii-1);
    end
end

b=zeros(K,M);       %Matrix für b
for jj=1:K
    for ii=1:M
        b(jj,ii)=1+2*r(jj+1,ii);
    end 
end

c=zeros(K,M);       %Matrix für c 
for jj=1:K
    for ii=1:M
        if ii==1
            c(jj,ii)=-2*r(jj+1,ii+1);	
        else
            c(jj,ii)=-r(jj+1,ii+1);
        end
    end
end 

d=zeros(1,M);         %Vektor für d
d1=zeros(1,M);
e=zeros(1,M);
U3=zeros(K+1,M+1);
U3(1,:)=uIC;          %IC
U3(:,M+1)=uR;         %BC

for jj=1:K
    for ii=1:M
        if jj==1
            d1(ii)=uIC(ii);
        else
            d1(ii)=d2(ii);
        end

        if ii==1
            e(ii)=2*r(jj+1,ii+1)*dx*cc(jj+1);         %Neumannrand
        elseif ii==M
            e(ii)=r(jj+1,ii+1)*U3(jj+1,ii+1);         %Dirichletrand
        else
            e(ii)=0;
        end   
        
        
        d(ii)=d1(ii)+e(ii)+dt*f(jj,ii);
    end
 
        %Sternumwandlung
        cistern(1)=c(jj,1)/b(jj,1);
        distern(1)=d(1)/b(jj,1);
    for ii=2:M
        cistern(ii)=c(jj,ii)/(b(jj,ii)-a(jj,ii)*cistern(ii-1));
        distern(ii)=(d(ii)-a(jj,ii)*distern(ii-1))/(b(jj,ii)-a(jj,ii)*cistern(ii-1));
    end
    %Rückwärtssubstitution
    u(M)=distern(M);
    for ii=M-1:-1:1
        u(ii)=distern(ii)-u(ii+1)*cistern(ii);
    end
        d2=u;
        U3(jj+1,1:M)=u;
end


%% Zeichnen

[XX,TT]=meshgrid(LR:dx:RR,t);
figure(1)
step=K/10;
surf_x=XX(1:step:end,1:1:end);
surf_t=TT(1:step:end,1:1:end);

if r<=0.5
subplot(2,2,1);
surf_u1=U1(1:step:end,1:1:end);
surf(surf_x,surf_t,surf_u1);
title('Wärmeleitung mit FTCS');
xlabel('Länge [dm]');               %Anfangsbedingung
ylabel('Zeit [s]');
zlabel('Temperatur [°C]');          %Temperaturverteilung
else 
subplot(2,2,1);
t1=text(0.0,0.6,'Für das FTCS-Verfahren sind die Stabilitätskriterien nicht erfüllt ','Color','r','FontWeight','bold','FontSize',12);
axis off
end 
   
    
subplot(2,2,3);
surf_u3=U3(1:step:end,1:1:end);
surf(surf_x,surf_t,surf_u3);
title('Wärmeleitung mit BTCS');
xlabel('Länge [dm]');               %Anfangsbedingung
ylabel('Zeit [s]');
zlabel('Temperatur [°C]');          %Temperaturverteilung

if Vari==1
       
    subplot(2,2,2);  
    text(0.0,0.2,['r = ', num2str(r(1,1))]);
    text(0.0,0.3,['a = ', num2str(aa(1,1))]);
    text(0.0,0.4,['K = ', num2str(K),';   dt = ', num2str(dt)]);
    text(0.0,0.5,['M = ', num2str(M),';   dx = ', num2str(dx)]);
    syms x t
    uIC=sin(pi*x);
    cc=1*t;
    uR=0.2*t;
    text(0.0,0.6,['IC: uIC(x) = ', char(uIC)]);
    text(0.0,0.7,['BC: u_x(0,t) = ', char(cc),';   u(1,t) = ', char(uR)]);
    text(0.0,0.8,sprintf('PDE: u_t=a*u_x_x'));
    t1=text(0.0,0.9,sprintf('Eingangswerte: '),'FontWeight','bold','FontSize',12);
    axis off

elseif Vari==2
    
    subplot(2,2,2);
    syms x t
    uIC=sin(pi*x);
    cc=1*t;
    uR=0.2*t;
    f=500*x+1*t;
    text(0.0,0.2,['r = ', num2str(r(1,1))]);
    text(0.0,0.3,['a = ', num2str(aa(1,1)),';   f = ', char(f)]);
    text(0.0,0.4,['K = ', num2str(K),';   dt = ', num2str(dt)]);
    text(0.0,0.5,['M = ', num2str(M),';   dx = ', num2str(dx)]);
    text(0.0,0.6,['IC: uIC(x) = ', char(uIC)]);
    text(0.0,0.7,['BC: u_x(0,t) = ', char(cc),';   u(1,t) = ', char(uR)]);
    text(0.0,0.8,sprintf('PDE: u_t=a*u_x_x'));
    t1=text(0.0,0.9,sprintf('Eingangswerte: '),'FontWeight','bold','FontSize',12);
    axis off
    
elseif Vari==3

    subplot(2,2,2);
    text(0.0,0.4,['K = ', num2str(K),';   dt = ', num2str(dt)]);
    text(0.0,0.5,['M = ', num2str(M),';   dx = ', num2str(dx)]);
    syms x t dt dx a
    uIC=sin(pi*x);
    cc=1*t;
    uR=0.2*t;
    f=500*x+1*t;   
    aa=0.05*t+0.05*x;
    r=a*dt/(dx*dx);
    text(0.0,0.2,['r = ', char(r)]);
    text(0.0,0.3,['a = ', char(aa),';   f = ', char(f)]);
    text(0.0,0.6,['IC: uIC(x) = ', char(uIC)]);
    text(0.0,0.7,['BC: u_x(0,t) = ', char(cc),';   u(1,t) = ', char(uR)]);
    text(0.0,0.8,sprintf('PDE: u_t=a*u_x_x'));
    t1=text(0.0,0.9,sprintf('Eingangswerte: '),'FontWeight','bold','FontSize',12);
    axis off
    
elseif Vari==4

    subplot(2,2,2);  
    text(0.0,0.2,['r = ', num2str(r(1,1))]);
    text(0.0,0.3,['a = ', num2str(aa(1,1))]);
    text(0.0,0.4,['K = ', num2str(K),';   dt = ', num2str(dt)]);
    text(0.0,0.5,['M = ', num2str(M),';   dx = ', num2str(dx)]);
    syms x t
    uIC=sin(pi*x);
    cc=1*t;
    uR=0.2*t;
    text(0.0,0.6,['IC: uIC(x) = ', char(uIC)]);
    text(0.0,0.7,['BC: u_x(0,t) = ', char(cc),';   u(1,t) = ', char(uR)]);
    text(0.0,0.8,sprintf('PDE: u_t=a*u_x_x'));
    t1=text(0.0,0.9,sprintf('Eingangswerte: '),'FontWeight','bold','FontSize',12);
    axis off

end
%{
subplot(2,3,1);
surf_u1=U1(1:step:end,1:1:end);
surf(surf_x,surf_t,surf_u1);
title('Wärmeleitung mit FTCS (inhomogen)');
xlabel('x-Achse');                  %Anfangsbedingung
ylabel('t-Achse');
zlabel('z-Achse');                  %Temperaturverteilung



subplot(2,3,3);
%surf(surf_x,surf_t,surf_u1-surf_u2);
title('Differenzen FTCS');
xlabel('x-Achse');                  %Anfangsbedingung
ylabel('t-Achse');
zlabel('z-Achse');                  %Temperaturverteilung



dddddddd=surf_u1-surf_u3;

subplot(2,3,5);
surf(surf_x,surf_t,surf_u1-surf_u3);
title('Differenz FTCS-BTCS');
xlabel('x-Achse');                  %Anfangsbedingung
ylabel('t-Achse');
zlabel('z-Achse');                  %Temperaturverteilung

%}