Solving Coupled system of differential equation with NDSolve

Question

I am solving a system with coupled system of differential equation to obtain the time evolution of two vectors. The problem stated here is equivalent to solving Landau Lifschitz Gilbert equationLLG Equation for two systems with a coupling between them. In my problem, m1[t] and m2[t] are two 3 component column vectors in cartesian coordinates. m1[t]={m1x,m1y,m1z} and m2[t]={m2x,m2y,m2z} respectively. m1[t] and m2[t] are acted upon by a term called effective field (HAnis).This HAnis is derived from an energy term

EAnis[r, \[Theta], \[CurlyPhi]] := (K1*(Sin[\[Theta]]^2) + 
           K2*(Sin[\[Theta]]^4) + K3*(Sin[\[Theta]]^6)*(Cos[6 \[CurlyPhi]]));
HAnis = Simplify[
   r*Grad[EAnis[r, \[Theta], \[CurlyPhi]], {r, \[Theta], \[CurlyPhi]},
      "Spherical"]]/MS;

Theta and phi are the polar and azimuthal angles in a spherical polar coordinate geometry. The time derivative of m1[t] and m2[t] is given by

cons1[t_] := -gamma*(10^(-9))*(-Hexch*Cross[m1[t], m2[t]] + 
     Cross[m1[t], (HAnis + Hext)]);
tGilbdamp1[t_] := alphag*Cross[m1[t], cons1[t]];
LLGS1 = {m1'[t] == cons1[t] + tGilbdamp1[t], 
   m1[0] == ToSphericalCoordinates[{1, 0, 0}]};

and

cons2[t_] := -gamma*(10^(-9))*(-Hexch*Cross[m2[t], m1[t]] + 
     Cross[m2[t], (HAnis + Hext)]);
tGilbdamp2[t_] := alphag*Cross[m2[t], cons2[t]];
LLGS2 = {m2'[t] == cons2[t] + tGilbdamp2[t], 
   m2[0] == ToSphericalCoordinates[{-1, 0, 0}]};

I tried to obtain numerical solution by

    sol1 = NDSolve[{LLGS1, LLGS2}, {m1, m2}, {t, tStart, tFinish}, 
   StartingStepSize -> tStep, 
   Method -> {"FixedStep", 
     Method -> "ExplicitRungeKutta"}, 
   AccuracyGoal -> 10, PrecisionGoal -> 100, MaxSteps -> \[Infinity]];

And I have stumbled upon the following error messages.

Cross::nonn1: The arguments are expected to be vectors of equal length, and the number of arguments is expected to be 1 less than their length.

Cross::nonn1: The arguments are expected to be vectors of equal length, and the number of arguments is expected to be 1 less than their length.

Cross::nonn1: The arguments are expected to be vectors of equal length, and the number of arguments is expected to be 1 less than their length.

General::stop: Further output of Cross::nonn1 will be suppressed during this calculation.

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.

Any help will be highly appreciated. The constants in the above equations are as follows:

gamma = 1.760859*(10^11);
alphag = (10^(-4));
tStep = (10^(-3));
tStart = 0;
tFinish = 8;
K1 = (3*(10^5));
K2 = (5*(10^5));
K3 = (5*(10^5));
Hext = {0, 0, 0};
Hexch = 1000;
MS = 238732.41;

Answer 1

If we correct all typos and add data for $\theta, \phi$, then the code works

gamma = 1.760859*10^11; 
alphag = 10^(-4); 
tStep = 10^(-3); 
tStart = 0; 
tFinish = 8; 
K1 = 3*10^5; 
K2 = 5*10^5; 
K3 = 5*10^5; 
Hext = {0, 0, 0}; 
Hexch = 1000; 
MS = 238732.41; EAnis[r_, \[Theta]_, \[CurlyPhi]_] := K1*Sin[\[Theta]]^2 + K2*Sin[\[Theta]]^4 + K3*Sin[\[Theta]]^6*Cos[6*\[CurlyPhi]]; 
HAnis = Simplify[r*Grad[EAnis[r, \[Theta], \[CurlyPhi]], {r, \[Theta], \[CurlyPhi]}, "Spherical"]]/MS; 
cons1[t_] := ((-gamma)*((-Hexch)*Cross[m1[t], m2[t]] + Cross[m1[t], HAnis + Hext]))/10^9; 
tGilbdamp1[t_] := alphag*Cross[m1[t], cons1[t]]; 
LLGS1 = {Derivative[1][m1][t] == cons1[t] + tGilbdamp1[t], m1[0] == ToSphericalCoordinates[{1, 0, 0}]}; 
cons2[t_] := ((-gamma)*((-Hexch)*Cross[m2[t], m1[t]] + Cross[m2[t], HAnis + Hext]))/10^9; 
tGilbdamp2[t_] := alphag*Cross[m2[t], cons2[t]]; 
LLGS2 = {Derivative[1][m2][t] == cons2[t] + tGilbdamp2[t], m2[0] == ToSphericalCoordinates[{-1, 0, 0}]}; 
sol1 = NDSolve[{LLGS1, 
    LLGS2} /. {\[Theta] -> 1, \[CurlyPhi] -> 0}, {m1, m2}, {t, tStart,
    tFinish}]

Plot[m1[t] /. sol1, {t, tStart, tFinish}]

Plot[m2[t] /. sol1, {t, tStart, tFinish}]