C 中 Runge-Kutta 四阶方案的二阶 ODE

Question

我正在使用以下代码求解一个二阶 ODE 的简谐运动。代码给了我一个正弦响应，但振幅越来越大，实际上应该始终保持不变。谁能指出我犯的任何错误？

#include<stdio.h>
#include<conio.h>

#define f(t,x1,x2) x2
#define g(t,x1,x2) -x1

int main()
{
  FILE *fp;

  float x1_0, x2_0, t0, tn, h, yn, k1, l1, k2, l2, k3, l3, k4, l4, k, l,
    x1_n, x2_n;
  int i, n;
  x1_0 = 3;
  x2_0 = 0;
  t0 = 0;
  tn = 100;
  h = 0.1;
  n = tn/h;


  /* Calculating step size (h) */
  printf("h = %f\n", h);
  /* Runge Kutta Method */
  printf("\nt0\tx1_0\tx1_n\tx2_n\tx2_n\n");
  for (i = 0; i < n; i++)
  {
    k1 = h * (f(t0, x1_0, x2_0));
    l1 = h * (g(t0, x1_0, x2_0));
    k2 = h * (f((t0 + h / 2), (x1_0 + k1 / 2), (x2_0 + l1 / 2)));
    l2 = h * (g((t0 + h / 2), (x1_0 + k1 / 2), (x2_0 + l1 / 2)));
    k3 = h * (f((t0 + h / 2), (x1_0 + k2 / 2), (x2_0 + l2 / 2)));
    l3 = h * (g((t0 + h / 2), (x1_0 + k2 / 2), (x2_0 + l2 / 2)));
    k4 = h * (f((t0 + h / 2), (x1_0 + k3 / 2), (x2_0 + l3 / 2)));
    l4 = h * (g((t0 + h / 2), (x1_0 + k3 / 2), (x2_0 + l3 / 2)));
    k = (k1 + 2 * k2 + 2 * k3 + k4) / 6;
    l = (l1 + 2 * l2 + 2 * l3 + l4) / 6;
    x1_n = x1_0 + k;
    x2_n = x2_0 + l;
    printf("%0.4f\t%0.4f\t%0.4f\t%0.4f\t%0.4f\n", t0, x1_0, x1_n, x2_0, x2_n);
    t0 = t0 + h;
    x1_0 = x1_n;
    x2_0 = x2_n;
  }

  /* Displaying result */
  printf("\nValue of t at x1 = %0.2f is %0.3f", tn, x1_n);

  return 0;
}

代码给了我一个正弦响应，但振幅越来越大，实际上应该始终保持不变。我试过更改“h”值和其他更改，但结果没有改变。谁能指出我犯的任何错误？