#include <cstdlib>
#include <iostream>
#include<cmath>
#include <fstream>

using namespace std;

int main(int argc, char *argv[])
{
    start:
    //a=algebrico
    //e=euler
    //rk=runge-Kutta
    ofstream myfile;
    ofstream myfile2;
    ofstream myfile3;
  
    myfile.open("eulerk.txt",ios::trunc);
    myfile2.open("eulerk2.txt",ios::trunc);
    myfile3.open("eulerk3.txt",ios::trunc);
    double x0,x,dx,i,n,e,ya,y,ye,L,yrk,p,xis,k1,k2,k3,k4;
    int fim;
    //declaracao da funcao f
    double f(double xis);
    double fa(double xis, double x0, double y);

    cout<<"\n Resolve a eq. diferencial dy = sqrt(x)dx\n\n";

    cout<<"insira o valor inicial x0\n";
    cin>>x0;
    cout<<"insira o valor inicial y(x0)\n";
    cin>>y;
    cout<<"insira o valor do incremento dx\n";
    cin>>dx;

//for(p=0;p<4;p++)
//{

    L=10.0;
    n=L/dx;
    ye = y;
    ya= y;
    yrk=y;
    x=x0;
    xis=x0;

    myfile<<x<<"\t";
    myfile<<ye<<"\n";

   myfile2<<x<<"\t";
   myfile2<<ya<<"\n";

   myfile3<<xis<<"\t";
   myfile3<<yrk<<"\n";

//método de euler
    for(i=0;i<n;i++){
                     ye = ye + f(x)*dx;
                     x=x+dx;
		     ya = fa(x,x0,y);
		       //f = cos(x);
                       //ya=sin(x);

                     cout<<x<<"\t";
                     cout<<ye<<"\t";
                     cout<<ya<<"\n";

                     myfile<<x<<"\t";
                     myfile<<ye<<"\n";
                     myfile2<<x<<"\t";
                     myfile2<<ya<<"\n";
    }

//metodo de runge-kutta
cout<<"\n";
cout<<"\n";
cout<<"\n";
cout<<"\n";
cout<<"\n";

 for(i=0;i<n;i++)
 {
k1=f(xis);
k2=f(xis+dx*1.0/2.0);
k3=f(xis+dx*1.0/2.0);
k4=f(xis+dx);
yrk=yrk+dx*(1.0/6.0)*(k1+2.0*k2+2.0*k3+k4);
xis=xis+dx;
cout<<xis<<"\t";
cout<<yrk<<"\n";
myfile3<<xis<<"\t";
myfile3<<yrk<<"\n";
}

//dx=dx/2.0;
//}

    cout<<"aperte 1 para terminar o programa ou qualquer tecla para voltar ao inicio\n";
    cin>>fim;

    if (fim==1){
       myfile.close();
       myfile2.close();
       myfile3.close();
       system("PAUSE");
       return EXIT_SUCCESS;
       }
    else
        goto start;
}



//definicao da funcao f para euler e rk4
double f(double x)
{
  double func;
  func = sqrt(x);
  return func;
}

//definicao da funcao fa solucao analitica
//da eq. y'= f(x)
double fa(double x, double x0, double y)
{
  double funca;
  funca = (2.0/3.0)*(sqrt(x)*x - sqrt(x0)*x0) + y;
  return funca;
}