Генератор эйлеровых гармоник

---

Привет,

У меня есть гармонический осциллятор, определяемый по :
f=-k/r + b/r^5 (центральная сила // r положение // и я буду называть v быстротой)
если я хорошо справился, то скорость будет такой: v_n+1 = v_n + Delta_t *( - k/m*r+ b/m*1/r^5)
Я называю v_nr = Delta_t *( - k/m*r+ b/m*1/r^5)

Итак, позиция r такова: r_n+1 = r_n + Delta_t * v_nr +1/2.(r_n+1 - r_n-1)

Мне нужно найти v(t) и r(t).
Сначала я попытаюсь решить эту проблему без члена b (который является членом возмущения)

Что я уже пробовал:

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#define M 3 // in kg
#define k 1.0  //  kg/sec^2
#define N 50   // interation

double T = 2*pi*(m/k)^0.5;

double main (void) {  // initiation
  double r=0.0; double rr=0.0;
  double v=10.0; double vr=0.0;
}

int main() {

    FILE *fichier;
    fichier = fopen("without b perturbation", "w+");
    int t=0 ;
    for(t = 0; t <= N; t++) {
        rr = r + T * v ;
        vr = v + T * k/m *r ;
        v = v + T * ((k * x / m) + (k * xint / m)) / 2.0;
        x = x + (h) * (v + vint) / 2.0;
	r = rr + T * (v+vr);
	v = v
	r=rr;
	fprintf(fichier, "%d\t%lf\n", t, yy);
	printf("%d\t%lf\n", t, yy);
    }
    fclose(fichier);
    return 0;
}

---