% Metodos Numericos 
% Ficha nº9 ex. 1
% Carlos Balsa - 2007/2008 

clear all
close all;

a=0; b=0.5;      % limites da integraçao 
h=0.1;           % passo
np=(b-a)/h+1;    % numero de pontos
t=[a:h:b]';      % pontos tk=t0+k*h
y0=1;            % condiçao inicial


% Metodo de Euler Simples
 y=zeros(np,1);   % iniciar vector soluçao 
 y(1)=y0;         % valor inicial

for k=1:np-1
        tk=t(k);
        yk=y(k);
        dy=fty(tk,yk);
        hdy=h*dy;
        y(k+1,1)=y(k,1)+hdy;      
end 
disp('Metodo de Euler'); 
y


% Metodo de Euler Modificado
z=zeros(np,1);   % iniciar vector soluçao 
z(1)=y0;         % Condiçao inicial

for k=1:np-1
    % previsao
        tk=t(k);
        yk=z(k);
        K1=fty(tk,yk);
        z(k+1)=z(k)+h*K1;
       
     % correcçao
        tk=t(k+1);
        yk=z(k+1);
        K2=fty(tk,yk);
        z(k+1)=z(k)+h*(K1+K2)/2;
end 
disp('Metodo de Euler Modificado'); 
z

% Metodo de Runge-Kutta de 4ª ordem
w=zeros(np,1);   % iniciar vector soluçao 
w(1)=y0;         % Condiçao inicial

for k=1:np-1
        tk=t(k);
        yk=w(k);
        k1=fty(tk,yk);
        
        tk=t(k)+h/2;
        yk=w(k)+h/2*k1;
        k2=fty(tk,yk);
        
        yk=w(k)+h/2*k2;
        k3=fty(tk,yk);
        
        tk=t(k)+h;
        yk=w(k)+h*k3;
        k4=fty(tk,yk);
        
        w(k+1)=w(k)+h/6*(k1+2*k2+2*k3+k4);
 end 
disp('Metodo de Euler Modificado'); 
w


% soluçao exacta
sol=zeros(np,1);
for k=1:np
    tk=t(k);
    sol(k)=sqrt(2-exp(-tk^2));
end

disp('Soluçao exacta'); 
sol

plot(t,y,'*:r'); hold on;
plot(t,z,'+-.g'); hold on;
plot(t,w,'o--c'); hold on;
plot(t,sol,'x-'); hold on;
grid on; hold on; 
xlabel('t_k');hold on; 
ylabel('y_k');hold on; 
legend('Euler','Euler Modificado','Runge-Kutta','Valores exactos',2); 
