Exemplo 2B - Sistema Predador Presa/Método RK4

%EXEMPLO 2B: SISTEMA PREDADOR PRESA/MÉTODO RK4
%IFPB 09/08/2017
%EXEMPLO DE APLICACAO - SISTEMA PREDADOR X PRESA
%X - TEMPO
%Y - PRESA
%Z - PREDADOR
clear all, close all, clc,
A=0; B=15; H=0.01; H1=H/2;
M=(B-A)/H;
X(1)=A; Y(1)=80; Z(1)=30;
for J=1:M
    x=X(J); y=Y(J); z=Z(J);
    K1=FP(x,y,z);
    L1=GP(x,y,z);
    K2=FP(x+H1,y+H1*K1,z+H1*L1);
    L2=GP(x+H1,y+H1*K1,z+H1*L1);
    K3=FP(x+H1,y+H1*K2,z+H1*L2);
    L3=GP(x+H1,y+H1*K2,z+H1*L2);
    K4=FP(x+H,y+H*K3,z+H*L3);
    L4=GP(x+H,y+H*K3,z+H*L3);
    Y(J+1)=Y(J)+(H/6)*(K1+2*K2+2*K3+K4);
    Z(J+1)=Z(J)+(H/6)*(L1+2*L2+2*L3+L4);
    X(J+1)=X(J)+H;
end
[X' Y' Z']
a=plot(X,Y,'k',X,Z,'b');set(a,'LineWidth',2);
legend('PRESA','PREDADOR',0);
xlabel('TEMPO')
pause, close
a=plot(Z(1),Y(1),'ro',Z,Y,'b'); set(a,'LineWidth',2);
legend('VALOR INICIAL','PREDADOR X PRESA',0);
xlabel('PREDADORES'), ylabel('PRESAS')