matlab - 实施 Radau-IIA 方法

Question

我正在寻找在 MATLAB 中执行 3 阶的 Radau-IIA 方法来解决 DAE 系统。这是我的系统（罗伯问题）：

y1'=-0.04*y1+10^4*y2*y3
y2'= 0.04*y1-10^4*y2*y3-3*10^7*y2^2
0=y1+y2+y3-1

初始值为

y1(0)=1 and y2(0)=y3(0)=0.

Runge-Kutta 画面看起来像这样

1/3  |  5/12  -1/12
1    |  3/4    1/4
_____|_______________
     |  3/4    1/4

（对不起，我不知道是否有可能在 StackOverflow 中编写 LaTeX 或其他东西）。

我不知道该怎么做..