功能: LobattoIII-B 5级8阶隐式龙格库塔法求解常微分方程

格式:
ODESolverByLobattoIIIB58(f, x0, y0, x, yName)
ODESolverByLobattoIIIB58(f, x0, y0, x, yName, xName)
ODESolverByLobattoIIIB58(f, x0, y0, x, yName, xName, eps)

f : 符号矩阵变量, 里面存储一阶导表达式
x0: 数值, 表示自变量初值
y0: 矩阵变量, 里面存储x0对应的y值
x : 矩阵变量, 表示待求解点位值
yName: 字符串变量, 里面存储多个变量名称,每个名称以都号分隔,其个数与y0一致
xName: x变量名称, 默认为空
eps  : 数值, 求解时控制的相对误差,默认为1e-15

说明:
1. 本函数主要求解一阶常微分初值问题.
2. f里表达式个数必须与y0里元素个数一致, 且一一对应.
3. yName里名称个数必须与y0里元素个数一致, 且一一对应.
4. x里面的值可以乱序, 返回的结果每行与x一一对应.
5. 本函数属于隐式求解算法，对于非刚性方程，建议使用显式求解算法。

参考: 
[1] Butcher J C. Numerical Methods for ordinary and differential equations[M]. John Wuley & Sons, 2016.

例子:
//已知x,y,t满足如下表达式, 现在需要求解t=[0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]时对应的x,y值
$$
\begin{cases}
x'&=-x
\\
y'&=-1000y
\\
x(0)&=1
\\
y(0)&=1
\end{cases}
$$
//这里执行如下代码
f = {SymMatrix<矩阵运算\SymMatrix>}([-x,
-1000*y]);
t0   = 0;
xy0 = [1,1];
t = [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1];
xy = ODESolverByLobattoIIIB58(f, t0, xy0, t, "x,y", "t")//回车得到如下结果
xy =
[ 0.90484048887440    2.6122763582E-26
  0.81873631030647    6.8239877719E-52
  0.74082576327693    1.7826141925E-77
  0.67032914581426    4.656680911E-103
  0.60654095200533    1.216453745E-128
  0.54882281153485    3.177713359E-154
  0.49659710109462    8.301065482E-180
  0.44934116372807    2.168467710E-205
  0.40658207825910    5.664636934E-231
  0.36789192645954    1.479759714E-256 ]

//可以将t带入下面理论表达式进行验证
$$
\begin{cases}
x &=\exp(-t)
\\
y &= \exp(-1000t)
\end{cases}
$$