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

float f1 ( float x, float y );
float f2 ( float x, float y );
void euler ( float f ( float x, float y ), float dx, int n, float x[], float y[] );
void ode_print ( int n, float x[], float y[] );

int main ( )
{
  float dx;
  int n;
  float x[51];
  float y[51];
//
//  Estimate 50 solution points along the curve defined by dy/dx = cos(x)
//  with starting point y(0.0) = 1.0, going from x = 0 to x = 5.
//
  dx = 0.1;
  n = 50; 
  x[0] = 0.0;
  y[0] = 1.0;

  euler ( f1, dx, n, x, y );

  printf ( "\n" );
  printf ( "Solution of dy/dx = cos(x)\n" );
  printf ( "\n" );

  ode_print ( n, x, y );
//
//  Now estimate 50 solution points along the curve defined by dy/dx = -20.0 * (x-1)*y
//  with starting point y(0.0) = 0.1, going from x = 0 to x = 2.
//
  dx = 0.04;
  n = 50; 
  x[0] = 0.0;
  y[0] = 0.1;

  euler ( f2, dx, n, x, y );

  printf ( "\n" );
  printf ( "Solution of dy/dx = -20*(x-1)*y\n" );
  printf ( "\n" );

  ode_print ( n, x, y );

  return 0;
}

float f1 ( float x, float y )
{
  float value;
  value = cos ( x );
  return value;
}

float f2 ( float x, float y )
{
  float value;
  value = -10.0 * ( x - 1.0 ) * y;
  return value;
}
 
void euler ( float f ( float x, float y ), float dx, int n, float x[], float y[] )
{
  float dydx;
  int i;

  for ( i = 0; i < n; i++ )
  {
    dydx = f ( x[i], y[i] );
    x[i+1] = x[i] + dx;
    y[i+1] = y[i] + dydx * dx;
  }
  return;
}

void ode_print ( int n, float x[], float y[] )
{
  int i;
  for ( i = 0; i <= n; i++ )
  {
    printf ( "%2d  %14g  %14g\n", i, x[i], y[i] );
  }
  return;
}