What methods can be used to numerically solve a system of linear non-differentiable equations?

I'm writing an implementation of Backwards Differentiation Formula, it's a multistep method of solving a system of ODEs. The four-step method looks like this:

$$y_{n+1} - y_n = h f(t_{n+1}, y_{n+1}) \\ y_{n+2} - \tfrac43 y_{n+1} + \tfrac13 y_n = \tfrac23 h f(t_{n+2}, y_{n+2}) \\ y_{n+3} - \tfrac{18}{11} y_{n+2} + \tfrac9{11} y_{n+1} - \tfrac2{11} y_n = \tfrac6{11} h f(t_{n+3}, y_{n+3}) \\ y_{n+4} - \tfrac{48}{25} y_{n+3} + \tfrac{36}{25} y_{n+2} - \tfrac{16}{25} y_{n+1} + \tfrac{3}{25} y_n = \tfrac{12}{25} h f(t_{n+4}, y_{n+4}) \\$$

At first I thought that I can solve it by quasi Newton-Raphson method, but it turns out that differential of $f$ can't be analytically found. Any suggestions?
I would appreciate pseudocode or code in any language if available.

-
Why not approximate the derivative by a secant? – copper.hat Jun 6 '12 at 5:29
Is there a generalization of secant method for the system of equations? – Andrew Jun 6 '12 at 7:59
Broyden's method (or some improvement) perhaps? – copper.hat Jun 6 '12 at 8:08
Thanks, checking it. Are you aware of any tutorials/code samples for it? – Andrew Jun 6 '12 at 8:49
For Broyden's method implemented in matlab, you can check phys.washington.edu/users/mforbes/projects/broyden – tibL Jun 6 '12 at 9:42