对于这样的ODEs难道你还想将它化为那个一阶显式微分方程组吗，别费劲了！
【解】
(1)好，我们同样的需要一组状态变量
(2).写出状态变量的一阶微分。只不过此时的x’’和y’’不再是显式的，我们使用fsolve进行数值求解
(3)编写Matlab代码
1. 2. 3. 4. 5. 6. 7. 8. 9.
[t,y]=ode45(@odefun,[0 3],[1 0 0 1]'); plot(t,y)
legend('x','x''','y','y''')
function dy=odefun(t,x)
fun=@(dxy)[dxy(1)*sin(x(4))+dxy(2)^2+2*x(1)*x(3)-x(1)*dxy(1)*x(4) x(1)*dxy(1)*dxy(2)+cos(dxy(2))-3*x(3)*x(2)]; options=optimset('display','off');
y=fsolve(fun,x([1,3]),options);%使用fsolve求解出x''和y''
10. dy=[x(2);y(1);x(4);y(2)];%状态变量一阶微分值 复制代码
================方法三================
方法二也有它的缺点，每一个点的x''和y''都需要独立计算，假如说数据量很大的话，那个叫难熬呀，并且有时可能是由于方程或者参数配置问题fsolve还有会出现罢工现象！
于是MathWorks为我们准备了ode15i函数，不过这个ode15i有点不好用，没有ode45等那么直接。
ode15i()调用格式如下
[y0mod,yp0mod,resnrm] = decic(odefun,t0,y0,fixed_y0,yp0,fixed_yp0, options...)
[T,Y,TE,YE,IE] = ode15i(odefun,tspan,y0mod,yp0mod,options...)
隐式ODEs不同于一般的显式ODEs，ode15i需要同时给出状态变量及其一阶导数的初值
(x0,dx0)，它们不能任意赋值，只能有n个是独立的，其余的需要隐式方程求解(使用decic函数)，否则出现矛盾的初值条件。
输入参数：
odefun：微分方程，注意此时的函数变量有t(自变量)，x(状态变量)，dx(状态变量一阶导数)等三个
t0：自变量初值，必须与tspan(1)相等
y0：状态变量初值，题目给出几个就写几个，没给的自己猜测一个
fixed_y0：与y0的尺寸一样，表示对应位置的初值是给定的还是猜测的，如果给定则1，否则0
yp0：状态变量一阶导数初值，其它同y0
fixed_yp0：同理fixed_y0，不过此时对应的是yp0了 options：其它选项，具体查看帮助
tspan：自变量的范围，其中tspan(1)必须与t0一致
输出参数：
大部分同理ode45，这里就不具体说明了
下面我们以实例说明，对于下面的IDE，求解x
【解】
(1)选取状态变量
(2)写出状态变量的一阶微分。只不过此时的x’’和y’’不再是显式的
(3)书写odefun，注意此时多了一个dx变量，就是x的一阶导数
1. 2. 3. 4.
odefun=@(t,x,dx)[dx(1)-x(2)
dx(2)*sin(x(4))+dx(4)^2+2*x(1)*x(3)-x(1)*dx(2)*x(4) dx(3)-x(4)
x(1)*dx(2)*dx(4)+cos(dx(4))-3*x(3)*x(2)];
(4)求解dx0
1. 2. 3. 4. 5. 6. 7. 8. 9.
t0=0;%自变量初值
%对于x0和dx0中的题目给出的初值，如实写，没有给出的任意写 x0=[1 0 0 1]';%本题初值x0的都给出了，所以必须是这个 %初值x0中题目给出的，对应位置填1，否则为0 fix_x0=ones(4,1);%本题中x0都给出了，故全为1
dx0=[0 0 1 1]';%本题中初值dx0一个都没有给出，那么全部任意写 %初值dx0中题目给出的，对应位置填1，否则为0 fix_dx0=zeros(4,1);%本题中dx0一个没有给出，故全部为0 [x02,dx02]=decic(odefun,t0,x0,fix_x0,dx0,fix_dx0);
(5)求解微分方程
1. solution=ode15i(odefun,[0 2],x02,dx02)%x02和dx02是由decic输出参数
虽然我们上面的分析是完全正确，但是好像Matlab就是罢工，郁闷死了
其实对一些简单的ODEs压根没有必要通过decic函数来求解未知的初值，我们可以直接人
工算出x02和dx02，根据下面的式子
我们演示下，根据题目x0=[x1,x2,x3,x4]=[1 0 0 1]，由第一个式子可以得出x1’=0，第三个式子有x3’=1，再联立第二四两个式子后有x2’=0，x4’=1，故dx0=[x1’ ,x2’,x3’,x4’]=[0 0 1 1]
看了下，隐式ODEs数值求解的确麻烦的些，但是和前面的ode45没有本质的区别，decic函数就是为了求解那些隐式中的状态变量的一阶导数，与我们方法二中的很相似，只不过这里直接调用罢了。
其实不管隐式还是显式ODEs，甚至DAE，都是可以使用ode15i求解，只不过此时将显式当隐式处理，或者说将代数约束看成一个隐式罢了！ 这一条等到我们在介绍DAE的时候再验证！
当然大部分人都不会傻到这么做，有现成的函数不用，何必自找麻烦呢！但是对于了解和学习，我们可以试试。
[教程] 微分代数方程(DAE)的Matlab解法
所谓微分代数方程，是指在微分方程中，某些变量满足某些代数方程的约束。假设微分方程的更一般形式可以写成
前面所介绍的ODEs数值解法主要针对能够转换为一阶常微分方程组的类型，故DAE就无法使用前面介绍的常微分方程解法直接求解，必须借助DAE的特殊解法。
其实对于我们使用Matlab求解DAE时，却没有太大的改变只需增加一个Mass参数即可。描述f(t,x)的方法和普通微分方程完全一致。
Matlab解微分方程(ODE+PDE)
