#include <stdio.h>
#include <math.h>
#include   <string.h>
#include   <stdlib.h>
#include   <alloc.h>
#define TRUE 1
#define FALSE 0
#define PI acos(-1)

int k_index;
float th=0.015;
float u_plant,y_plant,y_plant_old,ydot,ydot_old;
void rmodel();
float rmodel1(float);
void main (void)
{
  int i;

  y_plant=ydot=0.0;
  printf("t=0:0.015:7.5;\n");
  for(k_index=1; k_index<=501; k_index++)
   {
     u_plant=1.0;
     rmodel();
     printf("y(%d)=%f;\n",k_index,y_plant);
   }
  printf("\n plot(t,y);\n");
}

 void rmodel(void)
 {
     float x1;
     float g1,g2,g3,g4;

     x1 = y_plant;
     g1 = rmodel1(x1);
     g2 = rmodel1(x1+th/2*g1);
     g3 = rmodel1(x1+th/2*g2);
     g4 = rmodel1(x1+th*g3);

     y_plant += th/6*(g1+2*g2+2*g3+g4);
 }

 float rmodel1(float a)
 {
  float den;

  den = -1*a+0.5*a*a+u_plant;
  return(den);
 }