Numerische Integration einer DGL als C-Programm

#include <math.h>
#include <stdio.h>


float f(float x, float y);             // Rechte Seite der DGL  y'=f(x,y)


int main(void)
{
  const float x0 = 0;                  // Bereich, fuer den die Loesung bestimmt werden soll
  const float x1 = 0.7;
  const int   n  = 100;                // unterteile Bereich in n Schritte
  const float y0 = 1;                  // Startwert


  float dx = (x1-x0)/n;                // Schrittweite
  float dy;

  float y[n+1];                        // Tabelle fuer die zu berechnenden Funktionswerte


  y[0] = y0;                           // Der Funktionswert bei x0 ist bekannt


  for (int i=0 ; i<n ; i++)
  {
    dy = f(x0+i*dx, y[i]) * dx;

    y[i+1] = y[i] + dy;
  }


  for (int i=0 ; i<=n ; i++)           // Tabelle der Funktionswerte ausgeben
    printf("%f %f\n", x0+i*dx, y[i]);


  exit(0);
}



float f(float x, float y)
{
  float r;

  r = (1+y*y)*3*x;

  return r;
}

gcc -o dgl dgl.cc
dgl > data10  bzw. data50, data100, data500
gnuplot
> plot [0:0.8][0:10] tan(3*x**2/2+3.1415926/4), 'data500' w p, 'data100' w p, 'data50' w p, 'data10' w p