% Exercicio 1 p.15 do capitulo sobre EDPs
% Matematica Aplicada 2007/2008 - Carlos Balsa
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Resolver o exercicio proposto para a equaçao do calor usando
% o metodo das linhas.
% Depois de transformar a EDP num sistema de EDOs de primeira 
% ordem, aplicamos o metodo de Runge-Kutta de quarta ordem.
%
% Necessario ter o ficheiro dy.m no mesmo directorio,
% a funçao dy calcula o valor da derivada para os valores 
% de y e t requisitados no argumento da funçao
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
clear all;

c=1; % cosntante da equaçao do calor;

%% definiçao das linhas ao longo das quais queremos calcular a soluçao
dx=0.1;                 % passo espacial
a=0+dx;                 % ponto inicial 
b=1-dx;                 % ponto final
x=[a:dx:b]   ;          % valores de x 
nx=length(x);           % dimensao de x (nº de EDOs de 1ª ordem)

ti=0;            % tempo inicial
tf=1;            % tempo final
dt=dx^2/(2*c);   % passo de tempo
t=[ti:dt:tf] ;   % vector dos tempos
nt=length(t);    % dimensao do vector t


y=zeros(nt,nx);  % matriz das soluçoes iniciada a zero
y0=sin(x*pi);    % condiçoes iniciais
y(1,:)=y0;


% inicio da integraçao do sistema de EDOs 
% pelo metodo de Runge-Kutta de 4ª ordem
for k=1:nt-1

        tk=t(k);
        yk=y(k,:)';
        k1=dy(tk,yk,dx,nx,c); %  calculo da derivada atraves de funçao

        tk=t(k)+0.5*dt;
        yk=y(k,:)'+0.5*dt*k1;
        k2=dy(tk,yk,dx,nx,c);

        yk=y(k,:)'+0.5*dt*k2;
        k3=dy(tk,yk,dx,nx,c);

        tk=t(k)+dt;
        yk=y(k,:)'+dt*k3;
        k4=dy(tk,yk,dx,nx,c);
 
        media=dt/6*(k1+2*k2+2*k3+k4);

        y(k+1,:)=y(k,:)+media';
    end 
 

% grafico das linhas
 close all;
 waterfall(t,x,y'); hold on;  % grafico dos valores nas linhas de y
 view(45,35); hold on;        % ponto de vista
 colormap([0 0 0]);hold on;   % cor preta
 xlabel('t'); hold on; 
 ylabel('x'); hold on; 
 zlabel('u'); hold on; 
 
 