How can I pass a symbolic argument/function to Findroot or NDSolve and force them to not be evaluated until a numeric value is inputted?

Question

I have a 35X35 symbolic matrix that I would like to calculate the determinant of. With that being said, it takes way too long. As it is eventually going to be an equation in a system I feed to NDSolve, I was considering forcing it to be unevaluated until NDSolve takes it and "plugs" numbers in. Can someone tell me if this is possible with Mathematica?

My current implementation takes a system h(x)=0 and joins this system with $det\frac{\partial h}{\partial x}==0$ and attempts to solve using findroot. Call $det\frac{\partial h}{\partial x}=A$. Then with IC equal to an inital guess of solutions, I use

ndet[A_ /; MatrixQ[A, NumericQ]] := Det[N[A]];
sol = FindRoot[
  Join[h, {ndet[A]}] , IC]

As A is the determinant of a large (35x35) and symbolic matrix, I would like to only calculate A when the function Findroot is plugging in numbers.

I have a similar issue later on in the same code using the NDSolve, which I will omit here as the problem is similar.

As a side note, the findroot solution is used as an inital condition for a nonlinear differential-algabraic system later on.

Here is a workable snippet of the code up to FindRoot solution

(*enter parameters here*)
mug = 0.05;(*Jgenerator/Ms/R^2*)
Ms = 362.7;
Mus = 40;
u = Mus/Ms;(*mus/ms*)
V = 22.352;(*m/s driving speed*)
kt = 182087;(*tire stiffness*)
ks = 20053;(*suspension stiffness*)
omegac = 2*Pi*V*.005;
omega0 = Sqrt[ks/Ms];
omegat = Sqrt[kt/Ms];
vc = omegac/omega0;
v = omegat/omega0;
gr = .1617035985;(*road roughness*)
R = .03/2/Pi;
d = .01^2*7/omega0^2/R^2;(*2piGrV/omega0^2/R^2*)
mur = .1;(*maximize mur at .2 for compactness*)
cm = 148;(*mechanical damping in Ns/m*)
para = {\[Mu]r -> mur, xi -> cm/2/omega0/Ms, \[Mu]g -> mug};
(*\[Phi]0+\[Phi] expansion after M^-1. \
x[1]=\[Theta];x[2]=\[Phi];x[3]=\[CapitalPsi]s;x[4]=\[CapitalPsi]us*)
m[1, 1] = 1 + \[Mu]g + \[Mu]r*(1 + \[Eta]^2 + 2*\[Eta]*Cos[x[2]]);
m[1, 2] = \[Mu]r*\[Eta]*(\[Eta] + Cos[x[2]]) - \[Mu]g;
m[1, 3] = 1;
m[1, 4] = 0;
m[2, 1] = m[1, 2];
m[2, 2] = \[Mu]r*\[Eta]^2 + \[Mu]g;
m[2, 3] = 0;
m[2, 4] = 0;
m[3, 1] = 1;
m[3, 2] = 0;
m[3, 3] = 1 + u;
m[3, 4] = 0;
m[4, 1] = 0;
m[4, 2] = 0;
m[4, 3] = 0;
m[4, 4] = 0;
(mass,damping,stiffness matrix)
M = {{m[1, 1], m[1, 2], m[1, 3], m[1, 4]}, {m[2, 1], m[2, 2], m[2, 3],
     m[2, 4]}, {m[3, 1], m[3, 2], m[3, 3], m[3, 4]}, {m[4, 1], 
    m[4, 2], m[4, 3], m[4, 4]}};
Dm = {{2(xi + xie), -2 xie, 0, 0}, {-2 xie, 2xie, 0, 0}, {0, 0, 0, 
    0}, {0, 0, 0, 1}};
K = {{1, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, uv^2, -uv^2}, {0, 0, 0, vc}};
(external and nonlinear force)
Force = {[Mu]r[Eta]x[4](2x[3] + x[4])Sin[x[2]], -[Mu]r[Eta]*
    x[3]^2*Sin[x[2]], 0, n};
(substitute in harmonic functions)
EOM = (M . {x'[5], x'[6], x'[7], 0} + Dm . {x[5], x[6], x[7], x'[4]} +
    K . {x[1], x[2], x[3], x[4]} - Force)
sol = Solve[{EOM[[1]] == 0, EOM[[2]] == 0, EOM[[3]] == 0, 
   EOM[[4]] == 0}, {x'[4], x'[5], x'[6], x'[7]}]
EQ[1] = x'[5] /. sol[[1]];
EQ[2] = x'[6] /. sol[[1]];
EQ[3] = x'[7] /. sol[[1]];
EQ[4] = x'[4] /. sol[[1]];
eq[1] = (EQ[1] /. 
     x[2] -> [Phi]0) + (D[EQ[1], x[2]] /. x[2] -> [Phi]0)x[2];
eq[2] = (EQ[2] /. 
     x[2] -> [Phi]0) + (D[EQ[2], x[2]] /. x[2] -> [Phi]0)x[2];
eq[3] = (EQ[3] /. 
     x[2] -> [Phi]0) + (D[EQ[3], x[2]] /. x[2] -> [Phi]0)x[2];
eq[4] = (EQ[4] /. 
     x[2] -> [Phi]0) + (D[EQ[4], x[2]] /. x[2] -> [Phi]0)x[2];
(drift coefficients as a vector)
A = {x[5], x[6], x[7], eq[1], eq[2], eq[3], eq[4]} /. 
    n -> 0 /. -1 - [Mu]r + [Mu]rCos[[Phi]0]^2 -> -1 - [Mu]r
      Sin[[Phi]0]^2;
[CapitalSigma] = {0, 0, 0, Coefficient[eq[1], n]Sqrt[d], 
    Coefficient[eq[2], n]Sqrt[d], Coefficient[eq[3], n]Sqrt[d], 
    Coefficient[eq[4], n]
     Sqrt[d]} /. -1 - [Mu]r + [Mu]rCos[[Phi]0]^2 -> -1 - [Mu]r
      Sin[[Phi]0]^2;
Col = Outer[Times, [CapitalSigma], [CapitalSigma]];
numstates = 7;(number of state variable eqns)
(moment differential equations w/cumulant neglect)(define two 

differential operators here)
D1[f_, i_] := D[f, x[i]];
D2[f_, i_, j_] := D[f, {x[i], 1}, {x[j], 1}];
(g:moment generating function)
g[i_, j_, k_, n_, z_, y_, q_] = 
  x[7]^ix[6]^{j - i}x[5]^{k - j}x[4]^{n - k}x[3]^{z - n}*
   x[2]^{y - z}*x[1]^{q - y};
(generate a list of moments up to h-th order)
G = Flatten[
   Table[g[i, j, k, n, z, y, q], {q, 1, 2}, {y, 0, q}, {z, 0, y}, {n, 
     0, z}, {k, 0, n}, {j, 0, k}, {i, 0, j}]];
(define two Gaussian closure functions)
GClosure1[i_, j_, 
   k_] := [CapitalEpsilon][i][CapitalEpsilon][
     jk] + [CapitalEpsilon][j][CapitalEpsilon][
     ik] + [CapitalEpsilon][k][CapitalEpsilon][ij] - 
   2[CapitalEpsilon][i][CapitalEpsilon][j][CapitalEpsilon][k];
GClosure2[i_, j_, k_, 
   l_] := [CapitalEpsilon][ij][CapitalEpsilon][
     kl] + [CapitalEpsilon][i*k][CapitalEpsilon][
     jl] + [CapitalEpsilon][i*l][CapitalEpsilon][jk] - 
   2[CapitalEpsilon][i][CapitalEpsilon][j][CapitalEpsilon][
     k][CapitalEpsilon][l];
(close 3rd and 4th order moments with 1st and 2nd moments using 

Gaussian closure)
subs1 = Flatten[
   Table[x[i]x[j]x[k] -> GClosure1[x[i], x[j], x[k]], {i, 1, 7}, {j,
      1, 7}, {k, 1, 7}]];
subs2 = Flatten[
   Table[x[i]x[j]x[k]*x[l] -> GClosure2[x[i], x[j], x[k], x[l]], {i,
      1, 7}, {j, 1, 7}, {k, 1, 7}, {l, 1, 7}]];
(SUBs applies in order subs2,subs1,and Table[G[[i]][Rule]

[CapitalEpsilon][G[[i]]],{i,Length[G],1,-1}].The last substitutes 

[CapitalEpsilon][G[[i]]] for G[[i]],where [CapitalEpsilon][] 

denotes expected values)
SUBs = Join[subs2, subs1, 
   Table[G[[i]] -> [CapitalEpsilon][G[[i]]], {i, Length[G], 1, -1}]];
(moment equations.[CapitalEpsilon][G[[i]]][Rule]y[i][t] 

substitutes state variables y[i][t] for [CapitalEpsilon][G[[i]]])
SME = Flatten[
    Expand[Table[-Sum[D1[G[[i]], j]A[[j]], {j, 1, numstates}] - 
        1/2
         Sum[D2[G[[i]], j, k]*Col[[j, k]], {j, 1, numstates}, {k, 1, 
           numstates}], {i, 1, Length[G]}]] /. SUBs] /. 
   Table[[CapitalEpsilon][G[[i]]] -> y[i], {i, 1, Length[G]}];
Table[[CapitalEpsilon][G[[i]]] -> y[i], {i, 1, Length[G]}]
SME2 = SME /. Table[y[i] -> 0, {i, {1, 2, 3, 4, 5, 6, 7, 12}}];
y8sub = Flatten[Solve[{SME2[[11]] == 0}(/.[Phi]0->0), {y[8]}]] // 
  Simplify
y27sub = Flatten[Solve[{SME2[[27]] == 0}(/.[Phi]0->0), {y[27]}]] //
   Simplify
y10sub = Flatten[Solve[{SME2[[22]] == 0}(/.[Phi]0->0), {y[10]}]] //
   Simplify
y11sub = Flatten[Solve[SME2[[14]] == 0, y[11]]] // Simplify;
SMEsimplifieda = SME /. y8sub /. y11sub /. {y[13] -> -y[18]} /. Table[y[i] -> 0, {i, {1, 2, 3, 4, 5, 6, 7, 12}}] // Simplify; y29sub = Flatten[ Solve[SMEsimplifieda[[Length[SME`simplifieda]]] == 0, y[29]]]
SMEsimplifieda = SMEsimplifieda /. y29sub
(SME`simplifieda=SME/.y8sub/.y11sub/.y29sub/.{y[13]->-y[18]}/.Table[

y[i]->0,{i,{1,2,3,4,5,6,7,12}}]//Simplify;)
remaining = SparseArray[SME`simplifieda]["NonzeroPositions"] // Flatten
(SMEsimplified=SME/.y8sub/.y11sub/.y29sub/.{y[13]->-y[18]}/.Table[y[\ i]->0,{i,{1,2,3,4,5,6,7,12}}]; SMEsimplified=SME`simplified[[remaining]])
SMEsimplified = SMEsimplifieda[[remaining]]
Dimensions[%]
vars = Join[Table[y[i][[Eta]], {i, 9, 10}], 
  Table[y[i][[Eta]], {i, 14, 28}], 
  Table[y[i][[Eta]], {i, 30, Length[SME]}], {[Phi]0[[Eta]]}]
Length[vars]
varsd = D[vars, [Eta]];
subs = Join[
   Table[y[i] -> y[i][[Eta]], {i, 1, 
     Length[SME]}], {[Phi]0 -> [Phi]0[[Eta]]}];
AE = SME`simplified /. para /. subs;
AE2 = AE /. [Phi]0[[Eta]] -> 0
remainingAE = SparseArray[AE2]["NonzeroPositions"] // Flatten;
AE = AE[[remainingAE]]
Dimensions[%]
tt = Table[
     Coefficient[D[AE, [Eta]][[i]], varsd[[j]]], {i, 1, 
      Length[vars]}, {j, 1, Length[vars]}] /. 
    Join[Table[
      y[i][[Eta]] -> y[i], {i, 1, 
       35}], {[Phi]0[[Eta]] -> [Phi]0}] /. [Phi]0 -> 0 // Simplify
SparseArray[tt]["NonzeroPositions"]
SparseArray[tt]["AdjacencyLists"]
(tt=Grad[SME[[1;;Length[SME]]],Delete[Join[Table[y[i],{i,1,Length[

SME]}],{[Phi]0}],2]]/.Table[y[i][Rule]0,{i,{1,2,3,4,5,6,7,12}}]/.

Join[Table[y[i][[Eta]][Rule]y[i],{i,1,35}],{[Phi]0[[Eta]][Rule]

[Phi]0}]/.para/.[Phi]0[Rule]0//Simplify)
ndet[A_ /; MatrixQ[A, NumericQ]] := Det[N[A]];
(d`AE=Det[tt]//Simplify)
BF = AE /. 
     Join[Table[
       y[i][[Eta]] -> y[i], {i, 1, 
        35}], {[Phi]0[[Eta]] -> [Phi]0}] /. para /. [Phi]0 -> 0;
remainingBF = SparseArray[BF]["NonzeroPositions"] // Flatten;
BF = BF[[remainingBF]] // Simplify;
Eqn = Join[BF(,{d`AE})]
Length[Eqn]
vars2 = Join[Table[y[i], {i, 9, 10}], Table[y[i], {i, 14, 28}], 
  Table[y[i], {i, 30, Length[SME]}]]
Length[%]
IC = Join[
  Table[{vars2[[i]], .1}, {i, 1, 
    Length[vars2]}], {{xie, .1, .00000001, 10}}]
[Eta]0 = .9;
sol = FindRoot[
  Join[Eqn[[1 ;; Length[Eqn]]], {Det[
      tt]}] /. {[Eta] -> [Eta]0, [Phi]0 -> 0}, IC, 
  AccuracyGoal -> 14, PrecisionGoal -> 14]
Table[sol[[i]], {i, 1, Length[sol]}];
ICond = Join[% /. 
    Rule[x_, z_] -> x[0] == z, {[Eta][0] == [Eta]0}] /. 
  Table[y[i] -> ToExpression[StringJoin["y", ToString[i]]], {i, 1, 
    35}]

Edit; I have just edited my code snippet as there was a missing line