我正在使用下面的代码来解决一个简单的谐波运动，这是一个二阶常微分方程。代码给我一个正弦响应，但振幅越来越大，实际上应该保持不变。有人能指出我犯的任何错误吗？

#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”值和其他变化，但结果没有改变。有人能指出我犯的任何错误吗？This is how it should lookThis is how my code does it