f1[x_, y1_, y2_] := y2; 
      f2[x_, y1_, y2_] := -k/m*y1 (1 - l/(Sqrt[y1^2 + s^2]));
  
C2 = 1;  C1 = 1.  - C2; a2 = 1. /(2 C2); b2 = a2; 
  
h = 0.1; n = 100
  
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 = 1.45; y02 = 0; 
    RK1 = {y01}; RK2 = {y02}

s = 3; l = Sqrt[10. ]; m = 3. ; 
      k = 100

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

n

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 Exacta e Erro

Sol = 
      NDSolve[{y1'[t] == 
            y2[t], y2'[t] == -k/m*
              y1[t] (1 - l/(Sqrt[y1[t]^2 + s^2])), 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ÁFICO DA SOLUÇÃO EXACTA

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ÁFICO DA SOLUÇÃO APROXIMADA

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ÁFICO DO ERRO

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

Erro2

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

Show[Graf1, Graf2]

erro2par