@Joe Riel 

Hi Joe,

Thanks for your help.  I've ended up using 

for tt from 1 to max(t_steps) do
syseqs_ans(tt) := solve(eval(connecteqs union component_eqs,def_eqs)) union def_eqs; 
od:

Where connecteqs union component_eqs was my system of algebriac equations, and def_eqs contained my algebraic inputs (ie u1[1] = 3.5, u2[1] = 5.4, u1[2] = 3.3, u2[2] = 5.5 ....  .  The system solves very quickly (~100 systems/sec) with solve, but only about 2 systems/sec with fsolve.  I don't know why solve is so much faster than fsolve, but it works!

@Joe Riel 

Hi Joe.  Thanks for the tip.  I'll stop trying to do it this way then.

I have a system of (10+) algebraic equations (including soime 5th order polynomials) with about 3-5 inputs that I would like to simulate over 20000s (at 1s timestep).  I've tried a combination of map->fsolve(subs(control_inputs,EQ)), but solving this way takes too long.  I was hoping to use Simulate with AlgEquations, but you said I can't do this.  Do you have another method in mind that I could use?

Thanks Joe!