%function u=lampodiffE(r)
%% Lämpöyhtälö differenssimenetelmällä, Eksplsiittinen
%lampodiffE.m 
% Sauva olkoon 1:n pituinen
% Esim: u_t=u_xx, 
% alkuehto: u(x,0)=f(x)=sin(pi*x), 
% 0-reunaehdot.
clear; 
close all
n=6; h=1/(n-1)
x=0:h:1;
u0=sin(pi*x);   % Alkuarvot diskretiointipisteissä
%%  Muodostetaan iteraatiomatriisi M
r=0.25   % Varma valinta
%r=0.5;    % Suurin "sallittu"
%r=1;      % Aivan liian suuri
k=r*h^2
Tloppu=0.2;   % Loppuaika
t=0:k:Tloppu;
m=length(t);
u=zeros(n-2,m); % alueen sisällä n-2 x-jakopistettä
u(:,1)= u0(2:n-1)'; % jätetään päätyjen nollat pois 
M=lampomatr(n,r)
%%
%Iteroidaan u-matriisin sarakkeiksi eri aikatasot alkaen alkuarvoista u0, jotka
% jo yllä asetettiin sarakkeeksi 1.
for j=1:m-1
   u(:,j+1)= M*u(:,j);
end
 % Liitetään lopuksi matriisin ylle ja 
 % alle nollarivi, jotta saadaan reuna-arvotkin mukaan.
nollarivi=zeros(1,m);
u=[nollarivi;u;nollarivi];
%% Visualisointi
surf(x,t,u'); xlabel('x');ylabel('t'),zlabel('u')
title(['Lämpötilafunktion diffferenssiapproksimaatio',' r=',num2str(r)])
shg
figure
plot(x',u)
grid
xlabel('x');ylabel('u')
title(['Lämpötilan aikasiivujen kuvaajat hetkillä t=0:0.01:',num2str(Tloppu),' r=',num2str(r)])
shg
lampotilat=u'