Why is NDSolve not working for the given system?

Question

I'm trying to find the solution curves for a system of differential equations. In order to try to complete this task, I've attempted to use NDSolve, but I get an error that says that there are fewer dependent variables than equations, so the system is overdetermined. I've been working on this for hours and I can't seem to figure out how to fix the problem. I'm working on an SIR model with many variables so I apologize if it's difficult to read.

solcurve[\[ScriptCapitalS]m0_, \[CapitalEpsilon]m10_, \
\[CapitalIota]m10_, \[CapitalEpsilon]m20_, \[CapitalIota]m20_, \
\[ScriptCapitalS]h0_, \[CapitalEpsilon]h10_, \[CapitalIota]h1n0_, \
\[CapitalIota]h10_, Rh10_, \[CapitalEpsilon]h20_, \[CapitalIota]h2n0_, \
\[CapitalIota]h20_,Rh20_] := {\[ScriptCapitalS]m'[t] == \[Lambda]m - 
b*\[Beta]h1m*(\[ScriptCapitalS]m[t]*(\[CapitalIota]h1n[t] + \[CapitalIota]h1[
     t]))/\[CapitalNu]h - 
b*\[Beta]h2m*(\[ScriptCapitalS]m[
   t]*(\[CapitalIota]h2n[t] + \[CapitalIota]h2[
     t]))/\[CapitalNu]h - \[Mu]m*\[ScriptCapitalS]m[
  t], \[CapitalEpsilon]m1'[t] == 
b *\[Beta]h1m (\[ScriptCapitalS]m[
   t]*(\[CapitalIota]h1n[t] + \[CapitalIota]h1[
     t]))/\[CapitalNu]h - \[Gamma]m1 *\[CapitalEpsilon]m1[
  t] - \[Mu]m *\[CapitalEpsilon]m1[
  t] - \[Delta] *\[CapitalEpsilon]m1[t], \[CapitalIota]m1'[
t] == \[Gamma]m1*\[CapitalEpsilon]m1[t] - \[Mu]m*\[CapitalIota]m1[
  t] - \[Delta]*\[CapitalIota]m1[t], \[CapitalEpsilon]m2'[t] == 
b*\[Beta]h2m*(\[ScriptCapitalS]m[
   t]*(\[CapitalIota]h2n[t] + \[CapitalIota]h2[
     t]))/\[CapitalNu]h - \[Gamma]m2*\[CapitalEpsilon]m2[
  t] - \[Mu]m*\[CapitalEpsilon]m2[
  t] + \[Delta]*\[CapitalEpsilon]m1[t], \[CapitalIota]m2'[
t] == \[Gamma]m2*\[CapitalEpsilon]m2[t] - \[Mu]m*\[CapitalIota]m2[
  t] + \[Delta]*\[CapitalIota]m1[t], \[ScriptCapitalS]h'[
t] == \[Lambda]h - (
b*\[Beta]m1h*\[ScriptCapitalS]h[t]*\[CapitalIota]m1[
  t])/\[CapitalNu]h - (
b*\[Beta]m2h*\[ScriptCapitalS]h[t]*\[CapitalIota]m2[
  t])/\[CapitalNu]h - \[Mu]h*\[ScriptCapitalS]h[
  t], \[CapitalEpsilon]h1'[t] == (
b*\[Beta]m1h*\[ScriptCapitalS]h[t]*\[CapitalIota]m1[
  t])/\[CapitalNu]h - (\[Gamma]h + \[Mu]h)*\[CapitalEpsilon]h1[
  t], \[CapitalIota]h1n'[
t] == (1 - \[Phi])*\[Gamma]h*\[CapitalEpsilon]h2[
  t] - (q + \[Mu]h)*\[CapitalIota]h1n[t], \[CapitalIota]h1'[
t] == \[Phi]*\[Gamma]h*\[CapitalEpsilon]h1[
  t] - (q + \[Mu]h)*\[CapitalIota]h1[t], 
Rh1'[t] == 
 q*(\[CapitalIota]h2n[t] + \[CapitalIota]h2[t]) - \[Mu]h*
 Rh1[t], \[CapitalEpsilon]h2'[t] == (
b*\[Beta]m2h*\[ScriptCapitalS]h[t]*\[CapitalIota]m2[
  t])/\[CapitalNu]h - (\[Gamma]h + \[Mu]h)*\[CapitalEpsilon]h2 [
  t], \[CapitalIota]h2n'[
t] == (1 - \[Phi])*\[Gamma]h*\[CapitalEpsilon]h2[
  t] - (q + \[Mu]h)*\[CapitalIota]h2n[t], \[CapitalIota]h2'[
t] == \[Phi]*\[Gamma]h*\[CapitalEpsilon]h2[
  t] - (q + \[Mu]h)*\[CapitalIota]h2[t], 
  Rh2'[t] == 
   q*(\[CapitalIota]h2n[t] + \[CapitalIota]h2[t]) - \[Mu]h*
 Rh2[t], \[ScriptCapitalS]m[
0] == \[ScriptCapitalS]m0, \[CapitalEpsilon]m1[
0] == \[CapitalEpsilon]m10, \[CapitalIota]m1[
0] == \[CapitalIota]m10, \[CapitalEpsilon]m2[
0] == \[CapitalEpsilon]m20, \[CapitalIota]m2[
0] == \[CapitalIota]m20, \[ScriptCapitalS]h[
0] == \[ScriptCapitalS]h0, \[CapitalEpsilon]h1[
0] == \[CapitalEpsilon]h10, \[CapitalIota]h1n[
0] == \[CapitalIota]h1n0, \[CapitalIota]h1[
0] == \[CapitalIota]h10, 
  Rh1[0] == 
  Rh10, \[CapitalEpsilon]h2[
t] == \[CapitalEpsilon]h20, \[CapitalIota]h2n[
0] == \[CapitalIota]h2n0, \[CapitalIota]h2[
0] == \[CapitalIota]h20, Rh2[0] == Rh20}

Using solcurve above I then use NDSolve:

solucurve = 
 NDSolve[solcurve[1000000, 0, 10, 0, 15, 1500000, 0, 10, 90, 0, 0, 20,
10, 0], {\[ScriptCapitalS]m[t], \[CapitalEpsilon]m1[
t], \[CapitalIota]m1[t], \[CapitalEpsilon]m2[t], \[CapitalIota]m2[
t], \[ScriptCapitalS]h[t], \[CapitalEpsilon]h1[
t], \[CapitalIota]h1n[t], \[CapitalIota]h1[t], 
   Rh1[t], \[CapitalEpsilon]h2[t], \[CapitalIota]h2n[
    t], \[CapitalIota]h2[t], Rh20[t]}, {t, 0, 100}]

and these are the parameters: $(b=10;) (\text{$\beta $h1m}=0.8;) (\text{$\beta $h2m}=0.4;) (\text{$\beta $h2m}=0.6;) (\text{$\beta $m1h}=0.2;) (\text{$\beta $m2h}=0.7;) (\text{$\gamma $h}=0.4;) (\text{$\gamma $m1}=0.25;) (\text{$\gamma $m2}=0.5;) (\delta =3;) (\text{$\lambda $h}=0.7;) (\text{$\lambda $m}=2;) (\text{$\mu $h}=0.9;) (\text{$\mu $m}=0.1;) (\text{Nh}=100;) (q=0.3;) (y=1500000;) (\phi =0.4;)$

Answer 1

Anyway, let me do you a favor. It's easy to figure out where the typo is.

First one:

varInSystem = 
 Cases[solcurve[1000000, 0, 10, 0, 15, 1500000, 0, 10, 90, 0, 0, 20, 10, 0], _Symbol[t], 
   Infinity] // Union
(* {Rh1[t], 
 Rh2[t], \[ScriptCapitalS]h[t], \[ScriptCapitalS]m[t], Εh1[
  t], Εh2[t], Εm1[t], Εm2[
  t], Ιh1[t], Ιh1n[t], Ιh2[t], Ιh2n[
  t], Ιm1[t], Ιm2[t]} *)

var = {\[ScriptCapitalS]m[t], Εm1[t], Ιm1[
    t], Εm2[t], Ιm2[t], \[ScriptCapitalS]h[
    t], Εh1[t], Ιh1n[t], Ιh1[t], 
   Rh1[t], Εh2[t], Ιh2n[t], Ιh2[t], Rh20[t]};

Complement[var, varInSystem]
Complement[varInSystem, var]
(* {Rh20[t]} *)
(* {Rh2[t]} *)

Second one:

initInSystem = 
 Cases[solcurve[1000000, 0, 10, 0, 15, 1500000, 0, 10, 90, 0, 0, 20, 10, 0], _[0], 
   Infinity] // Union
(* {Rh1[0], 
 Rh2[0], \[ScriptCapitalS]h[0], \[ScriptCapitalS]m[0], Εh1[
  0], Εm1[0], Εm2[0], Ιh1[
  0], Ιh1n[0], Ιh2[0], Ιh2n[0], Ιm1[
  0], Ιm2[0]} *)

init = var /. t -> 0
(* {\[ScriptCapitalS]m[0], Εm1[0], Ιm1[
  0], Εm2[0], Ιm2[0], \[ScriptCapitalS]h[
  0], Εh1[0], Ιh1n[0], Ιh1[0], 
 Rh1[0], Εh2[0], Ιh2n[0], Ιh2[0], Rh2[0]} *)

Complement[init, initInSystem]
(* {Εh2[0]} *)