Just to clear something up about my previous post...
My matrix is actually 217x217 or 200x200, not 47089x47089 or 40000x40000. So I do not think my matrix size is the problem. What I think is the problem is that fsolve is taking my 40000 or 47089 elements and squaring that to make a jacobian... Does that make sense?
I am still trying to come up with a way to produce the Jacobian (yes, mathematically it's straightforward, but expressing it in a numerical way cleanly is another question.)
-John