2nd Order ODE by Runge-Kutta 4th order scheme in C

Question

I am using the following code to solve a simple harmonic motion which is a 2nd Order ODE. The code gives me a sinusoidal response, but the amplitude keeps getting bigger and bigger which in actual should stay constant throughout. Can anyone point out any mistake that I have made?

#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;
}

The code gives me a sinusoidal response, but the amplitude keeps getting bigger and bigger which in actual should stay constant throughout. I have tried changing the "h" values and other changes, but the result is not changing. Can anyone point out any mistake that I have made? This is how it should look This is how my code does it

Answer 1

There may be other issues, but one is in the way k4 and l4 are calculated:

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)));
//                  ^^^              ^^^              ^^^

Those values shouldn't be halved¹.

---

^{1) See e.g. https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods and the modified code here: https://godbolt.org/z/WKcdo4cav}