The numerical solution of a (very ugly) set of integro-diferential equations

I've been having fun playing around with a simple model for reaction in a co-flow fluidized bed (I'm afraid I'm a chemical engineer, not a computer scientist). Unfortunately, simple physical assumptions often lead to terrifying equations. I've reduced the problem to solving the following set of simultaneous integro-differential equations for the 2 unknown functions $g(s)$ (on $s\in[s_0,s_{max}]$) and $C(z)$ (on $z\in [0,L]$), where $f(C,s)$ is a known function and $\alpha$ and $\tau$ are constants:

$$-\frac{dC}{dz} = \alpha \int_{s_0}^{s_{max}} f(C(z),s) g(s) ds$$ $$g(s) = \frac{\exp\left(-\frac{L}{\tau}\int_{s_0}^{s} \frac{ds}{\int_{0}^{L}f(C(z),s)dz}\right)}{\frac{\tau}{L} \int_{0}^{L}f(C(z),s)dz}$$ and the one BC: $$C|_{z = 0} = C_0$$

So far I've tried the brute-force approach: finite differences. I divided the intervals $[0,L]$ and $[s_0,s_{max}]$ into $n$ points, and defined the $2n$ unknown variables $C_1, C_2,\ldots,C_n$ and $g_1, g_2,\ldots,g_n$, which are just the values of my unknown functions at these points. I then approximately the derivatives and integrals above using FD expressions, and attempted to solve the system of $2n$ non-linear algebraic equations using $\texttt{fsolve}$ in MATLAB.

Perhaps predictably, MATLAB doesn't like this, and fsolve struggles to give a solution for medium-sized $n$ ($n=100$ returns 'no solution found', and $n = 250$ is still running after half an hour.) My questions would be:

• Are there any ways of solving large systems of simultaneous non-linear equations, that are more efficient than $\texttt{fsolve}$. Preferably I'd be looking for an algorithm already implemented in MATLAB.
• Even more optimistically, are there any other possible approaches to solving these types of highly non-linear integro-differential equations?

Thanks

• I may be misunderstanding things, but how is g an unknown function, given you have an (integral) formula in terms of C? – origimbo Feb 26 '16 at 13:42
• @origimbo Yes, you're 100% right, in reality I only have one unknown function, C (z), as I could substitute the second equation directly into the first and eliminate g(s). I didn't realise this until after I wrote this post, because from a physical perspective its intuitive to think of C (a concentration profile) and g (a solids concentration distribution) as the two unknown functions in the system. The aforementioned substition would lead to one single very ugly id equation, which I would still have no idea how to solve. :) – tom Feb 26 '16 at 14:25