% Metodo de previsao correcçao de 2ª ordem
% Exercici 1 (b) da ficha nº 6
% Matematica Aplicada 2008/2009
% Carlos Balsa

clear all

t0=0; tf=2; h=0.25; 
np=(tf-t0)/h+1;
y=zeros(np,1);
t=[t0:h:tf]';
y0=1; y(1)=y0;
y(2)=0.94427; % valor de y1

for k=2:np-1
    % previsao
        tkm1=t(k-1); tk=t(k);  
        ykm1=y(k-1); yk=y(k);  
        dykm1=-2*tkm1*ykm1^2; dyk=-2*tk*yk^2; 
        y(k+1)=y(k)+h/2*(3*dyk-dykm1);
  
       
     % correcçao
        dykm1=dyk; 
        tk=t(k+1);
        yk=y(k+1);
        dyk=-2*tk*yk^2; 
        y(k+1)=y(k)+h/2*(dyk+dykm1);
end 

 plot(t,y,'r+:')
y
