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.