#Get the tight coupling equations for the scalar modes

pb(t):=0;
drag(t):= drag_t;
opac(t):=opc(t)/eps;
photbar(t):=photbar_t;
rhob(t):=4*rhog(t)/3/R(t);
polter(t):=polter_t;

eqs:={simplify(diff(pig(t),t)-dpig),simplify(diff(E2(t),t)-dE2),simplify(diff(B2(t),t)-dB2),simplify(diff(J_3(t),t)-dJ_3)};


makeseries:=proc(v,x,ord)
 local Res,i;
 Res:=0;
 for i from 0 to ord do
   Res:=Res + cat('s',v)[i](t)*x^i;
 end do;
 v(t)=series(Res,x,ord); 
end;

solvevars := {Delta,pig,J_3,B2,B3,E2,E3,J_3,J_4};

series_subs:=map(makeseries,solvevars,eps,4);

subs(series_subs,eqs);

eqs:=map(series,%,eps,3);
map(coeff,eqs,eps,-1);
spig[0](t):=0;
sDelta[0](t):=0;
sE2[0](t):=0;
sB2[0](t):=0;
sB3[0](t):=0;
sB3[1](t):=0;
sJ_3[0](t):=0;
sE3[0](t):=0;
sJ_3[0](t):=0;
sJ_3[1](t):=0;
sJ_4[0](t):=0;

#sJ_3[1](t):=0;
#sE3[1](t):=0;
#sJ_4[1](t):=0;
#sJ_4[2](t):=0;


eqs0:=map(coeff,eqs,eps,0);

solve(eqs0,{spig[1](t),sE2[1](t),sB2[1](t)});