(*Método de Runge - 
          Kutta para sistemas de duas equações*)

f1[x_, y1_, y2_] := y1 (a - b y2); 
      f2[x_, y1_, y2_] := -y2 (c - d  y1);
  
C2 = 0.5;  C1 = 1 - C2; a2 = 1/(2 C2); b2 = a2; 
  
h = 0.05; n = 200
  
V1[x_, y1_, y2_] := 
        f1[x + a2 h, y1 + b2 h f1[x, y1, y2], 
          y2 + b2 h f2[x, y1, y2]] ; 
      V2[x_, y1_, y2_] := 
        f2[x + a2 h, y1 + b2 h f1[x, y1, y2], 
          y2 + b2 h f2[x, y1, y2]] ;
  
x0 = 0; 
    y01 = 20; y02 = 5; RK1 = {y01}; RK2 = {y02}

a = 1. ; b = 0.1; c = 1.0; 
      d = 0.1;

(* Cálculo da solução aproximada*)

Do[xi = x0 + i *h;  
      ynovo1 = RK1[[-1]] + 
          h (C1  f1[xi, RK1[[-1]], RK2[[-1]]] + 
                C2 V1[xi, RK1[[-1]], RK2[[-1]]]); 
      ynovo2 = RK2[[-1]] + 
          h (C1  f2[xi, RK1[[-1]], RK2[[-1]]] + 
                C2 V2[xi, RK1[[-1]], 
                    RK2[[-1]]]); 
AppendTo[RK1, 
        ynovo1]; AppendTo[RK2, ynovo2], {i, 0, n - 1}]

Solução  pelo Comando NDSolve

Sol = 
      NDSolve[{y1'[t] == y1[t] (a - b y2[t]), y2'[t] == 
            y2[t] (-c + d  y1[t]), y1[0] == y01, 
          y2[0] == y02}, {y1, y2}, {t, 0, n*h}]

Exacta1 = Table[y1[i*h] /. Sol[[1]], {i, 0, n}]; 
    Exacta2 = Table[y2[i*h] /. Sol[[1]], {i, 0, n}]; 
    Erro1 = Exacta1 - RK1; Erro2 = Exacta2 - RK2;

GRÁFICOS DA SOLUÇÃO OBTIDA PELO NDSOLVE

Graf = 
      Plot[Evaluate[{y1[t], y2[t]} /. Sol], {t, 0, n*h}, 
        PlotStyle -> {{AbsoluteThickness[0.05], 
              RGBColor[0, 0, 1]}, {AbsoluteThickness[0.05], 
              RGBColor[1, 0, 0]}}]

GRÁFICOS DA SOLUÇÃO OBTIDA PELO MÉTODO DE RUNGE - KUTTA

RK1par = Table[{(i - 1)*h, RK1[[i]]}, {i, 1, n}]; 
    Graf1 = ListPlot[RK1par, 
        PlotStyle -> {RGBColor[0, 0, 1]}]

RK2par = Table[{(i - 1)*h, RK2[[i]]}, {i, 1, n}]; 
    Graf2 = ListPlot[RK2par, 
        PlotStyle -> {RGBColor[1, 0, 0]}]

Show[Graf1, Graf2]

Show[Graf, Graf1, Graf2]

GRÁFICOS DO ERRO

erro1par = Table[{(i - 1)*h, Erro1[[i]]}, {i, 1, n}]; 
    Graf1 = ListPlot[erro1par, 
        PlotStyle -> {RGBColor[0, 0, 1]}]

erro2par = Table[{(i - 1)*h, Erro2[[i]]}, {i, 1, n}]; 
    Graf2 = ListPlot[erro2par, 
        PlotStyle -> {RGBColor[1, 0, 0]}]

Show[Graf1, Graf2]