Falkner skan solution with different beta value
 

I am trying to write falkner skan code with varying value of beta(2*m/m+1), where m tell us about velocity profile. My code looks like...

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

using namespace std;

double function_f1(double z,double z1,double z2,double eta);
double function_f2(double z,double z1,double z2 ,double eta);
double function_f3( double z,double z1,double z2,double eta,double beta);

int main(){

double beta,deta,res=10.0;
double *z,*z1,*z2,*z2new;
double *eta;
int e,max_e=7001,iter;
double k[4],m[4],n[4];
double A,A_improved;


eta=new double [max_e+1];
z=new double [max_e+1];
z1=new double [max_e+1];
z2=new double [max_e+1];
z2new=new double [max_e+1];

//initial values
eta[0]=0.0L;
deta= 0.001L;
z[0]=0.0L;
z1[0]=0.0L;
z2[0]=0.322L;

A=z2[0];
double *diff,*err,*q;
diff=new double [50];
err =new double [50];
q =new double [50];
beta = -0.0 ;
iter=0;

for (int iter=0;iter<12;iter++)
{
for (e=0;e<max_e;e++)
{
k[0]= function_f1(z[e],z1[e],z2[e],eta[e]);
m[0] = function_f2(z[e],z1[e], z2[e],eta[e]);
n[0] = function_f3(z[e],z1[e],z2[e],eta[e], beta);

k[1]= function_f1(z[e]+0.5*k[0]*deta,z1[e]+0.5*m[0]*deta,z2[e]+0.5*n[0]*deta,eta[e]+0.50*deta);
m[1] = function_f2(z[e]+0.5*k[0]*deta,z1[e]+0.5*m[0]*deta,z2[e]+0.5*n[0]*deta,eta[e]+0.50*deta);
n[1] = function_f3(z[e]+0.5*k[0]*deta,z1[e]+0.5*m[0]*deta,z2[e]+0.5*n[0]*deta,eta[e]+0.50*deta,beta);

k[2]= function_f1(z[e]+0.5*k[1]*deta,z1[e]+0.5*m[1]*deta,z2[e]+0.5*n[1]*deta,eta[e]+0.50*deta);
m[2] = function_f2(z[e]+0.5*k[1]*deta,z1[e]+0.5*m[1]*deta,z2[e]+0.5*n[1]*deta,eta[e]+0.50*deta);
n[2] = function_f3(z[e]+0.5*k[1]*deta,z1[e]+0.5*m[1]*deta,z2[e]+0.5*n[1]*deta,eta[e]+0.50*deta,beta);

k[3]= function_f1(z[e]+k[2]*deta,z1[e]+m[2]*deta,z2[e]+n[2]*deta,eta[e]+0.50*deta);
m[3] = function_f2(z[e]+k[2]*deta,z1[e]+m[2]*deta,z2[e]+n[2]*deta,eta[e]+0.50*deta);
n[3] = function_f3(z[e]+k[2]*deta,z1[e]+m[2]*deta,z2[e]+n[2]*deta,eta[e]+0.50*deta,beta);

z[e+1] = z[e]+ (deta/6.0)* ( k[0] + 2.0*k[1] +2.0*k[2] + k[3] ) ;
z1[e+1] = z1[e]+ (deta/6.0)* ( m[0] + 2.0*m[1] +2.0*m[2] + m[3] );
z2[e+1] = z2[e]+ (deta/6.0)* ( n[0] + 2.0*n[1] +2.0*n[2] + n[3] );

eta[e+1] =eta[e]+ deta;

//diff[iter]= z1[1088]; //random point near 1.0 choosen from the curve

}

ofstream outfile;
outfile.open("falknerbeta0.0.dat");

for (e=0;e<max_e;e++)
{
cout<< eta[e] <<" "<<z[e]<<" "<<z1[e]<<" "<<z2[e]<<endl;
outfile<< eta[e] <<" "<<z[e]<<" "<<z1[e]<<" "<<z2[e]<<endl;
}
outfile.close();

/*//newton raphson
q[iter]= A;
err[iter] = (1.0L - z1[1088]);

if (iter == 1) {A=A+ 0.05;}

else

if (iter == 3)
if (err[2] <err[1])
A_improved=q[2]-((q[3]-q[2])*err[2]/(err[3]-err[2]));
else
{A_improved=q[1]-((q[3]-q[1])*err[1]/(err[3]-err[1]));}

else
{ A_improved=q[iter-1]-((q[iter]-q[iter-1])*err[iter-1]/(err[iter]-err[iter-1]));}


{q[iter]=q[iter-1]-((q[iter]-q[iter-1])*err[iter-1]/(err[iter]-err[iter-1]));}

if (fabs(err[iter]> 0.0001)){ z2[0]= q[iter];}
else {exit(0);}*/
z2[0]=z2[0]+0.001;

}

delete z2new;
delete z2;
delete z1;
delete z;

return 0;

}



double function_f1(double z,double z1,double z2,double eta){

double f1;
f1= z1;
return f1;
}


double function_f2(double z,double z1,double z2,double eta){

double f2;
f2= z2;
return f2;
}


double function_f3( double z,double z1,double z2,double eta,double beta){

double f3;
f3= -0.5*z*z2 - beta * (1.0L-z1*z1);
return f3;
}


i am using RK4 with values try to refine using newton raphson...

My question is
1.if i am doing it correctly. Is there any better way to write this code in c++.
2. for blasius profile it is giving correct result but become unstable at higher beta value(beta =1.0, 0.5 ..) .why?(i dont know)
3. the commented part is newton raphson ...i tried to refine my rk4 solution with NR but could not implement it correctly...any one please try

thanks !