I am trying to implement an implicit Runge-Kutta method for
$$f(t,y)=-4t^2y^2+2y =y'$$
Do the $k_l$ in the description have to be unique ? In my example I got to solve a quadratic equation to solve for $k_1$ which would give me two values. Any help would be appreciated.
You have to chose the solution that is closer to $y_n$, so that $y_{n+1}-y_n=O(h)$. Usually there should only be one such solution that is close to the Euler step $y_n+hf(t_n,y_n)$.
For instance for the implicit Euler method (as it is the simplest one, not for its numerical properties), you have to solve $$ z=y_n+hf(t_{n+1},z)=y_n+h(−4t_{n+1}^2z^2+2z)\\ z=y+h(-4t^2z^2+2z)\\ (4ht^2)z^2+(1-2h)z+(-y)=0\\ $$ where the solution formulas give \begin{aligned} z&=\frac{-(1-2h)\pm\sqrt{(1-2h)^2+16ht^2y}}{8ht^2}\\[.8em] &=\frac{2y}{(1-2h)\pm\sqrt{(1-2h)^2+16ht^2y}} \end{aligned} Now one solution is $$ z=\frac{-(1-2h)-\sqrt{(1-2h)^2+16ht^2y}}{8ht^2} $$ and for $h<0.5$ this is of size $O(1/h)$. For the other solution use the second solution formula to find \begin{align} z&=\frac{2y}{(1-2h)+\sqrt{(1-2h)^2+16ht^2y}}\\ &=\frac{y}{1-2h+4ht^2y+O(h^2)} \end{align} which is indeed close to $y=y_n$ by $O(h)$. With the first term of the geometric series one finds further that $z=y(1+2h-4ht^2y+O(h^2))=y+hf(t,y)+O(h^2)$, which is close to the forward Euler step.
