% Varmeligning (sid 925 i Kreuzig)
clc;
clear all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Parametrar
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
L           =     1;    % Lengde er lik 1
c           =     1;    % Konstant c   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nods        =   6;             % Antall noder (i ROM)
h           =   L/(Nods-1);    % Intervallbredd (ROM)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N           =   20;     % Antall tidsskritt
k           =   h^2;     % Skrittlengde i tid, slik at r=1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r           =   k/(h^2); %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

UU=zeros(Nods,N+1);      % Minnesallokering;
                         % rader = Noder ROM, Kolloner = Tidsskritt

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Matriser
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A=zeros((Nods-2),(Nods-2));  % Matrise A
B=zeros((Nods-2),(Nods-2));  % Matrise B
% slik at UU(2:(Nods-1),j+1)=inv(B)*(A * UU(2:(Nods-1),j) + Randcerdier)
A=diag((2-2*r)*ones(Nods-2,1)) + diag(r*ones(Nods-3,1),1)...
    + diag(r*ones(Nods-3,1),-1);
B=diag((2+2*r)*ones(Nods-2,1)) - diag(r*ones(Nods-3,1),1)...
    - diag(r*ones(Nods-3,1),-1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Initial og randverdier
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
UU(1,1:(N+1))       =   zeros(1,N+1);    % Randverdi x=0
UU(Nods,1:(N+1))    =   zeros(1,N+1);    % Randverdi x=1
%
XX      =   [0:h:L]';        % 0 <= x <= 1
UU(:,1) =   sin(pi*XX);      % Initialverdi

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Numerisk lösning
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j=1:N
    Randverd=zeros((Nods-2),1);
    Randverd(1,1)=r*( UU(1,j) + UU(1,j+1) );
    Randverd((Nods-2),1)=r*( UU(Nods,j) + UU(Nods,j+1) );
    %
    UU(2:(Nods-1),j+1)=inv(B)*...
        ( A * UU(2:(Nods-1),j) + Randverd ); 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Plot 1:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1) % grid
TTid=[0:k:k*N]';
mesh(TTid,XX,UU)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Plot 1:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Opner en avi-fil med navnet Eksempel_film.avi
figure(2)
mov = avifile('Eksempel_film2.avi','compression','Cinepak','fps',15);

for i=1:(N+1)
    
    % Plotter den numeriske lösningen
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    hold off % Sletter forige figuren
    % plotte med röd (r) strek (-) å sirkel (o) '-or'
    plot(XX,UU(:,i),'r') 

    % Plotter den eksakte lösningen 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    hold on % Sparer forige figur
    x=[0:(L/100):L]; % x-verdier 
    ui_exakt=sin(pi*x)*exp(-pi^2*((i-1)*k));
    plot(x,ui_exakt,'-.','LineWidth',2)
    
    % Axis setter skalerer akslarna i figuren 
    axis ([0 L 0 1]); 
    
    %%%
    F = getframe(gca); 
    mov = addframe(mov,F);  % Legg til figur til film
    %%%
    
    pause(0.02)
    
end
mov = close(mov);  % Steng film

% Slut