% Runge-Kutta de 4ª ordem
 
format short;
 
clear all
close all
a=0; 
b=10;
h=0.1;
n=3;
np=(b-a)/h+1;
c1=1; c2=10;
y=zeros(np,n);
load rk_ex2ci;
y=temp(1:4,:);


t=[a:h:b]';
%y0=[1 1 1];
%y(1,:)=y0;
c2=10;



for k=4:np-1
   %Predictor
       tk=t(k);
        yk=y(k,:)'; % y_(k)
        dy3=[(-c1*yk(1)) (c1*yk(1)-c2*yk(2)) (c2*yk(2))]';
   
        tk=t(k-1);
        yk=y(k-1,:)'; % y_(k-1)
        dy2=[(-c1*yk(1)) (c1*yk(1)-c2*yk(2)) (c2*yk(2))]';

        tk=t(k-2);
        yk=y(k-2,:);
        dy1=[(-c1*yk(1)) (c1*yk(1)-c2*yk(2)) (c2*yk(2))]';

        tk=t(k-3);
        yk=y(k-3,:);
        dy0=[(-c1*yk(1)) (c1*yk(1)-c2*yk(2)) (c2*yk(2))]';

        media=h/24*(55*dy3-59*dy2+37*dy1-9*dy0);

        y(k+1,:)=y(k,:)+media';
       
        
      %Corrector
         tk=t(k+1);
         yk=y(k+1,:);
         dy4=[(-c1*yk(1)) (c1*yk(1)-c2*yk(2)) (c2*yk(2))]';
         media=h/24*(9*dy4+19*dy3-5*dy2+dy1);
         y(k+1,:)=y(k,:)+media';
     end 
 
 
% grafico
plot(t,y(:,1),'r*:'); hold on; 
plot(t,y(:,2),'b+:'); hold on;
plot(t,y(:,3),'k.-'); hold on;
xlabel('x_k'); hold on; 
ylabel('t_k'); hold on; 
 legend('y_1','y_2','y_3'); hold on; 
 

