Algorithm 196: Muller's method for finding roots of an abitrary function
Robert D. Rodman · Communications of the ACM · 1963
ZERSOL finds the simple zeros of the solution Y1 (Y0) of the set of m first order differential equations Yj = Fj(YO, Y1, ... , Ym). h is the step of integration, epsi the error with which the zeros are to be determined (assuming no error in the process of integration).F(YS, j, v) is a procedure which calculates the functions Fj, taking the arguments from the array YS and leaving the results in v.The search for zeros stops when Y0 > f.The zeros are stored as elements of the array Z. MR is a 4 X 4 matrix with the coefficients of a Runge-Kutta method.For example MR may be row-wise 0.5, 1, 0.5, 0, 1 -a, 1 --a, 1 --a, 0.5, l+a, 1 +a, 1 +a, 0,~,$, 0.5, 0.5, where a = sqrt(2); begin real v, r, d; integer ], s, n, ]6; array Q[l:m], YS[0:m], YAL [0:m], YT[i:m], MR[l:4,1:4]; switch S := NOZ, ZER; n:=l; for d := h while YI[0] ~ f do begin s := 1; Ri: forj := 1 step 1 until m do begin Q[j] := 0.0; YS[j] := Y/[j]; YT