(*  This notebook solves the semiclassical backreaction equations using Method 1.  C0 = N*λ and for N identical scalar fields the stress-tensor for the fields is proportional to it.  

I'll use A instead of the B I had Robert use.  A = B/3.  *)

(* Nf is the number of conformally invariant scalar fields. *)

αq = Nf/(2880*Pi^2)

βq = αq

(* The order reduced equations have these quantities in them.  There are no derivatives of hte scale factor on the right hand side of the equations.  *)

Abar = A*(1+(16*Pi*λ/3)*(2*αq+βq))

Bbar = (1+8*Pi*λ*βq/3)/3

Cbar = 8*Pi*λ*A^2*βq

(* Note that to bring this into the current notation, Nf λ = C0 *)

Abar = Abar/.{Nf*λ->C0}

Bbar = Bbar/.{Nf*λ->C0}

Cbar = Cbar/.{Nf*λ->C0}

(* This is the equation for α'^2.  *)

al1p2 = Abar/al^2-1+ al^2*Bbar+Cbar/al^6

al2p = -Abar/α[τ]^3 + Bbar*α[τ] + Cbar/α[τ]^7

(* This gives alpha at the bounce*)

s1 = Solve[(al1p2/.{A->1/2,C0->10^2})==0,al]

albounce = N[Part[s1,4,1,2]]

(* Check that α' is zero at this value of α *)

al1p2/.{al->albounce,A->1/2,C0->10^2}

(* Check that we have a minimum rather than a maximum. *)

al2p/.{α[τ]->albounce,A->1/2,C0->10^2}

(* Now solve the equation for α''.  *)

Eqmethod1 = NDSolve[{α''[τ]==(al2p/.{A->1/2,C0->10^2}), α'[0]==0,α[0] ==albounce},α,{τ,-10,10}]

Plot[α[t]/.Eqmethod1,{t,-10,10},PlotRange->All]

Save["C:\\Users\\anderson\\Research\\cosmology\\Order Reduction\\Conformally_invariant_fields\\Plots_for_the_paper\\Solutions_Method_1.m",Eqmethod1]