% Exercicio Proposto no capitulo sobre EDPs
% Matematica Aplicada 2007/2008 - Carlos Balsa
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Resolver a equaçao do calor usando um 
% esquema de diferenças finitas implicito.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
clear all;

c=1;
% variaveis espaciais
xa=0; xb=1; dx=0.25;
x=[xa:dx:xb];
nx=length(x)-2; % nº pontos a resolver no espaço

% variaveis temporais
ti=0; tf=0.3; dt=0.1;
t=[ti:dt:tf];
nt=length(t)-1; % nº de ptos a resolver no tempo

% condiçoes iniciais em t0
u0=sin(pi*x);

% condiçoes de fronteira 
 ua=zeros(nt,1); % em a
 ub=zeros(nt,1); % em b
 
% matriz com todos os pontos da malha
U=zeros(nt+1,nx+2);
U(1,:)= u0;
U(2:nt+1,1)=ua;
U(2:nt+1,nx+2)=ub;
 
gama=c*dt/dx^2; % constante em cada equaçao
N=nx*nt; % numero total de pontos a resolver

%% Geraçao da matriz dos coeficientes
d0=-(1+2*gama)*ones(nx,1) ;            
d1=gama*ones(nx-1,1) ; 
A=diag(d0)+diag(d1,1)+diag(d1,-1); 

%% Geraçao dos termos independente dos sistemas
b=zeros(N,1)
for j=1:nx;
    b(j)=-u0(j); % condiçoes iniciais
end

l=0;
for j=1:nx:N-nx+1;
    l=l+1;
    b(j)=b(j)-gama*ua(l); % condiçoes de fronteira em xa
end

l=0;
for i=nx:nx:N;
    l=l+1;
    b(j)=b(j)-gama*ub(l); % condiçoes de fronteira em xb
end

% resoluçao dos varios sistemas
for k=0:nt-1
    jj=k*nx
    y=A\b(jj+1:jj+nx);
    U(2+k,2:nx+1)=y'; 
    if (k~=nt-1)
    b(jj+nx+1:jj+2*nx)=b(jj+nx+1:jj+2*nx)-y; %actualiza b
    end
end

% grafico da superficie suluçao
 close all;
 mesh(t,x,U'); hold on;  % superficie 3d a partir dos pontos
 view(35,25);  hold on;  % ponto de vista
 xlabel('t');  hold on; 
 ylabel('x');  hold on; 
 zlabel('u');  hold on; 