Needs["LinearAlgebra`MatrixManipulation`"]
Relaxation[A_List, B_List, u_List, SlowFast_List,
           SlowVars_List, FastVars_List, 
           BCs_List, IV_List] := Block[
	{ns, nf, StateVars, T, dT, IVs, RHS, z0, z,
	 Ai, Bi, ABi, OneSolu, Solu, dti, i, j, m},
	
	m = Length[A];
	ns = Length[SlowFast[[1]]];
	nf = Length[SlowFast[[2]]];
	StateVars = Table[If[MemberQ[SlowFast[[1]], i],
	                     SlowVars[[Position[SlowFast[[1]], 
	                                        i][[1, 1]]]],
	                     FastVars[[Position[SlowFast[[2]], 
	                                        i][[1, 1]]]]
	                    ], {i, 1, ns + nf}];

	T = BCs;
	dT = {};
	IVs = {IV};
	Solu = {};
	z0 = Array[z, {nf}]; 
	For[i=1, i<=m, i++,
		ABi = ReduceMatrix[A[[i]], B[[i]], SlowFast];
		Ai = ABi[[1]];
		Bi = ABi[[2]];
	
		RHS = Table[Part[Ai . StateVars +
                                 Evaluate[Bi . u /. {t -> t + T[[i]]}],
		                 SlowFast[[2, j]]],
		            {j, 1, nf}];

		OneSolu = DSolve[Join[Table[D[FastVars, t][[j]] == 
		                            RHS[[j]], {j, 1, nf}],
		                      Table[FastVars[[j]] ==
		                            z0[[j]],
		                            {j, 1, nf}] /. {t -> 0}],
		                 FastVars, t] /. Table[z0[[j]] -> 
		                                       IVs[[i, j]],
                                               {j, 1, nf}];

	   OneSolu = Simplify[ComplexExpand[OneSolu]];

	   Solu = Join[Solu, {Table[OneSolu[[1, j, 2]],
	                            {j, 1, nf}]} ];

	   dti = If[!VectorQ[BCs[[i+1]]], BCs[[i+1]] - T[[i]],
                    If[Length[BCs[[i+1]]]==1,
                       BCs[[i+1, 1]],
                       Part[Solve[Evaluate[BCs[[i+1, 1]] == 
                                           BCs[[i+1, 2]] /. OneSolu[[1]]],
                                  t], 1, 1, 2]
                   ]
               ];

      If[MemberQ[dti, -1]&&MemberQ[dti, ArcSin[_]],
         dti = -Delete[dti, Part[Position[dti, ArcSin[_]],
                                 1, 1]] Pi - dti];
      
      T[[i+1]] = T[[i]] + dti;              
      
      dT = Join[dT, {dti}];

	  IVs = Join[IVs, {Simplify[Solu[[i]] /. {t -> dT[[i]]}]}]; 
	               
    ];
    
    {Solu, dT, T, IVs}
]
ReduceMatrix[A_List, B_List, SlowFast_List] := Block[
	{ABSlow, ABFast, A1, B1, 
	 NullVecA, NullVecB, i},

	ABSlow = Part[Transpose[Join[Transpose[A], 
	                             Transpose[B]]], 
	              SlowFast[[1]]];
	ABFast = Part[Transpose[Join[Transpose[A], 
	                             Transpose[B]]], 
	              SlowFast[[2]]];
	NullVecA = Table[0, {Length[A[[1]]]}];
	NullVecB = Table[0, {Length[B[[1]]]}];
	
	A1 = A;
	B1 = B;
	
	For[i=1, i<=Length[ABFast], i++,
		If[MemberQ[ABSlow, ABFast[[i]]],
		   A1[[SlowFast[[2, i]]]] = NullVecA];
		If[MemberQ[ABSlow, ABFast[[i]]],
		   B1[[SlowFast[[2, i]]]] = NullVecB];
		];
	
	{A1, B1}
]
Decoupling[FastSolus_List, A_List, B_List, SlowFast_List,
           SlowVars_List, u_List, T_List] := Block[
    {Xi, F, fi, i, j},

	F = {};
    For[i=1, i<=Length[A], i++,
    	Xi = Table[If[MemberQ[SlowFast[[1]], j],
	                  SlowVars[[Position[SlowFast[[1]], 
	                                     j][[1, 1]]]],
	                  FastSolus[[1, i,
	                             Position[SlowFast[[2]], 
	                                      j][[1, 1]]]]
	                  ], 
	              {j, 1, Length[A[[i]]]}];
		
		fi = {Simplify[Part[A[[i]] . Xi +
                               Evaluate[B[[i]] . u /. {t -> t + T[[i]]}],
		                    SlowFast[[1]]] ],
		      FastSolus[[2, i]]};
		F = Append[F, fi];
		];
	F
]
Averaging[SlowEqus_List, Ts_] := Block[{AveF, i, j},
		
	AveF = Sum[Integrate[SlowEqus[[i, 1]],
                             {t, 0, SlowEqus[[i, 2]]}],
		       {i, 1, Length[SlowEqus]}];
	
	Simplify[AveF/Ts]
]