NotebookAutoSave->True;

use Jeff Freidberg's function for <σv>
also use physical expression for (Subscript[τ, α]), taken from Betti ad Freidberg, PoF B 4, 1465 (1992)

note that we have defined:
(Subscript[τ, eq])=c1 τEi;
	(Subscript[τ, Ee])= c2 τEi;
	(Subscript[τ, α loss])=c4 τEi.

Eα=3.5 (10^3); (* [keV] *)

(* from Freidberg *)
(* n in (m^(-3)), T in keV *)
kB=0.052;
eV=1.6022 (10^(-19));
cB = ((10 )^(-1))eV kB;

σvJF[T_?NumericQ]=Exp[k0σ+k1σ Log[T]+k2σ ((Log[T])^2)+k3σ ((Log[T])^3)+k4σ ((Log[T])^4)];
k0σ=-60.4593;
k1σ=6.1371;
k2σ=-0.8609;
k3σ=0.0356;
k4σ=-0.0045;

σvfun[x_?NumericQ]=σvJF[x];


τα[Te_?NumericQ,n_?NumericQ]:=1.17 (10^18)((Te^(3/2))/n);
τei[Te_?NumericQ,Ti_?NumericQ,n_?NumericQ]:=(((78 ((n (10^(-20)))/(((((Te+Ti)/2)))^(3/2)))))^(-1));

(* We assign input profiles for density and temperatures, so that we can study the effect of the profile on ignition. Profiles are kept fixed in time, i.e. they are picked once and for all at the beginning of each simulation. Ion and electron temperature profiles do not need to be the same.
To avoid numerical issues, functions do not go to 0 at the edge, but to a finite value tedge or nedge. *)

tedge = (10^(-5));
nedge = (10^(-5));

tempμν[r_?NumericQ,μ_?NumericQ,ν_?NumericQ]:=(1-tedge)(((1-(r^ν)))^μ)+tedge (* this is Th *)
nηθ[r_?NumericQ,η_?NumericQ,θ_?NumericQ]:=(1-nedge)(((1-(r^θ)))^η)+nedge(* this is nh *)

(* Something needs to be made up for nα. We can use an arbitrary profile, or try and make one up that's derived from the system equations. We first try getting the nα profile from the α evolution equation. Observe that we solve for nαhat, we still take time-derivatives of nα[t] *)

(* Define an intermediate function and use funny variable names just in case; lval stands for local value. *)
fnα[r_?NumericQ,ηη_?NumericQ,θθ_?NumericQ,μμi_?NumericQ,ννi_?NumericQ,μμe_?NumericQ,ννe_?NumericQ,nlval_?NumericQ,Tilval_?NumericQ,Telval_?NumericQ,τEαlval_?NumericQ]:=(1/4)(((nlval nηθ[r,ηη,θθ]))^2)σvfun[Tilval tempμν[r,μμi,ννi]]/((1/(τα[Telval tempμν[r,μμe,ννe],nlval nηθ[r,ηη,θθ]]))+(1/(τEαlval)));


nαh[r_?NumericQ,ηη_?NumericQ,θθ_?NumericQ,μμi_?NumericQ,ννi_?NumericQ,μμe_?NumericQ,ννe_?NumericQ,nlval_?NumericQ,Tilval_?NumericQ,Telval_?NumericQ,τEαlval_?NumericQ]:=fnα[r,ηη,θθ,μμi,ννi,μμe,ννe,nlval,Tilval,Telval,τEαlval]/fnα[0,ηη,θθ,μμi,ννi,μμe,ννe,nlval,Tilval,Telval,τEαlval]

Notice that integrals are over the area, so in r dr dθ (here θ is the angular variable). For convenience and since nothing depends on the angle, we’ll do an integral in r dr only.
This also means that we are neglecting the effect of curvature (can put those back in later, if so desired).

(* For reference *)
(* H-mode-like profiles *)

μH=0.5;
νH=1.5;

ηH=0.25;
θH=1.25;

(* L-mode-like profiles *)

μL=1.5;
νL=2.5;

ηL=1.25;
θL=2.25;

Now we need to define the RHS of each differential equation (after integration in r dr). It seems reasonable to divide the nα equation by nαhat before integrating. This leaves only one messy integral.
Contrary to the original SF file, we leave the profile parameters in the definitions, so that if necessary one can iterate over parameters.

f1b[μ_?NumericQ,ν_?NumericQ,η_?NumericQ,θ_?NumericQ]:= NIntegrate[nηθ[r,η,θ]tempμν[r,μ,ν]r,{r,0,1}];
f2b[μ_?NumericQ,ν_?NumericQ,η_?NumericQ,θ_?NumericQ]:= NIntegrate[((nηθ[r,η,θ])^2)((tempμν[r,μ,ν])^(1/2))r,{r,0,1}]
f3b[μ_?NumericQ,ν_?NumericQ,η_?NumericQ,θ_?NumericQ]:= NIntegrate[(1/(τα[tempμν[r,μ,ν],nηθ[r,η,θ]]))r,{r,0,1}]

(* Original 0-D equations, for reference *)
rhs1[n_,nα_,Ti_,Te_,c2_,τEi_,Shi_,She_]:=Shi-1.5 ((n Ti)/(τEi))+1.5 ((n(Te-Ti))/(τei[Te,Ti,n]))
rhs2[n_,nα_,Ti_,Te_,c2_,τEi_,Shi_,She_]:=She-1.5 ((n Te)/(c2 τEi))+1.5 ((n(Ti-Te))/(τei[Te,Ti,n]))+((nα)/(τα[Te,n] ))Eα-cB (n^2) (Te^(1/2))
rhs3[n_,nα_,Ti_,Te_,c2_,c4_,τEi_,Shi_,She_]:=((n^2)/4)σvfun[Ti]-((nα)/(τα[Te,n]))-((nα)/(c4 τEi ))

(* complete 1-D equations *)
(* ion energy equation *)
rhs1b[n_?NumericQ,nα_?NumericQ,Ti_?NumericQ,Te_?NumericQ,c2_?NumericQ,τEi_?NumericQ,Shi_?NumericQ,She_?NumericQ,η_?NumericQ,θ_?NumericQ,μi_?NumericQ,νi_?NumericQ,μe_?NumericQ,νe_?NumericQ]:= Shi-1.5 ((n Ti)/(τEi))f1b[μi,νi,η,θ]+1.5 NIntegrate[ FullSimplify[((n nηθ[r,η,θ](Te tempμν[r,μe,νe]-Ti tempμν[r,μi,νi]))/(τei[Te tempμν[r,μe,νe],Ti  tempμν[r,μi,νi],n nηθ[r,η,θ]]))r],{r,0,1}]

(* electron energy equation *)
rhs2b[n_?NumericQ,nα_?NumericQ,Ti_?NumericQ,Te_?NumericQ,c2_?NumericQ,c4_?NumericQ,τEi_?NumericQ,Shi_?NumericQ,She_?NumericQ,η_?NumericQ,θ_?NumericQ,μi_?NumericQ,νi_?NumericQ,μe_?NumericQ,νe_?NumericQ]:=She-1.5 ((n Te)/(c2 τEi))f1b[μe,νe,η,θ]-1.5 NIntegrate[ FullSimplify[((n nηθ[r,η,θ](Te tempμν[r,μe,νe]-Ti tempμν[r,μi,νi]))/(τei[Te tempμν[r,μe,νe],Ti  tempμν[r,μi,νi],n nηθ[r,η,θ]]))r],{r,0,1}]+nα NIntegrate[(1/(τα[Te  tempμν[r,μe,νe],n nηθ[r,η,θ]] ))nαh[r,η,θ,μi,νi,μe,νe,n,Ti,Te,c4 τEi]r,{r,0,1}]Eα-cB (n^2) (Te^(1/2)) f2b[μe,νe,η,θ]

(* α density equation - remember: this is divided by nαh *)
rhs3b[n_?NumericQ,nα_?NumericQ,Ti_?NumericQ,Te_?NumericQ,c4_?NumericQ,τEi_?NumericQ,η_?NumericQ,θ_?NumericQ,μi_?NumericQ,νi_?NumericQ,μe_?NumericQ,νe_?NumericQ]:=((n^2)/4) NIntegrate[((((nηθ[r,η,θ])^2)σvfun[Ti tempμν[r,μi,νi]])/(nαh[r,η,θ,μi,νi,μe,νe,n,Ti,Te,c4 τEi]))r,{r,0,1}]-nα  (n/(Te^(3/2))) f3b[μe,νe,η,θ]-((nα)/(c4 τEi ))

(* set up the profile parameters *)

(* for example, use L-mode profiles *)

(* temperature profiles first *)

μitest=μL;
νitest=νL;

μetest=μL;
νetest=νL;

(* then density profiles *)

ηtest = ηL;
θtest = θL;

(* ratios between energy confinement times ans central density *)
c2loc=1;
c4loc=1;
nloc= (10^20);

(* ion energy confinement time *)
τEiloc=4;

(* electron and ion heating powers *)
Sheloc=0;
Shiloc=2. (10^20)



eqdiff1NLfull:=-1.5f1b[μitest,νitest,ηtest,θtest]nloc  Ti'[t]+rhs1b[nloc,nα[t],Ti[t],Te[t],c2loc,τEiloc,Shiloc,Sheloc,ηtest,θtest,μitest,νitest,μetest,νetest];

eqdiff2NLfull:=-1.5f1b[μetest,νetest,ηtest,θtest] nloc  Te'[t]+rhs2b[nloc,nα[t],Ti[t],Te[t],c2loc,c4loc,τEiloc,Shiloc,Sheloc,ηtest,θtest,μitest,νitest,μetest,νetest];

eqdiff3NLfull:=-nα'[t]+rhs3b[nloc,nα[t],Ti[t],Te[t],c4loc,τEiloc,ηtest,θtest,μitest,νitest,μetest,νetest];


tmax=30 τEiloc;

solNL=NDSolve[{eqdiff1NLfull==0,eqdiff2NLfull==0,eqdiff3NLfull==0,Ti[0]==0.002,Te[0]==0.001,nα[0]==0},{Ti[t],Te[t],nα[t]},{t,0,tmax}];
Ti1NL[t_]=First[Chop[Ti[t]/.solNL]];
Te1NL[t_]=First[Chop[Te[t]/.solNL]];
nα1NL[t_]=First[Chop[nα[t]/.solNL]];

tt=Te1NL["Domain"];
tdomain=Last[Last[tt]];

Tifin=Re[Ti1NL[tdomain]]
Tefin=Re[Te1NL[tdomain]]
nαfin=Re[nα1NL[tdomain]]

Plot[Ti1NL[t],{t,0,tmax},PlotRange->All]
Plot[Ti1NL[t],{t,0,τEiloc},PlotRange->All]