Format long;
N=(b-a)/h;
y=zeros(N+1,1);
y(1)=y0;
x=a:h:b;
var=findsym(f);
for i=2:N+1
   fx=Funval(f,var(1),x(i));
   gx=y(i-1)+h*fx-varvec(2);
   y(i)=NewtonRoot(gx,-10,10,eps);%用牛顿法得出下步的迭代值
end
format short;

yt[n+1] = y[n] + h * f(x[n], y[n]);
y[n+1] = y[n] + h * f(x[n+1], yt[n+1]);

3楼: Originally posted by sudo at 2012-02-23 13:29:11:
恩，如果不要求通用的方法，只需要针对特定函数的话，也可以先化简迭代式，或者可能可以直接导出y[n+1]的公式什么的~

4楼: Originally posted by 微尘、梦想 at 2012-02-23 13:38:29:
那个MATLAB程序就是用牛顿法得到的下一个值，不过，不知道牛顿法的算法是什么呀，我只知道用牛顿法解方程的公式：
x[n+1]=x[n]-f(x[n])/f`(x[n])
这是一种切线法，从一端向一个方向逼近方程的根，可是跟这个联 ...

   fx=Funval(f,var(1),x(i)); %里面的x(i)就相当于x[n+1]，var(1)是一个符号变量（其实是等下要求解的y[n+1]），把它们代入到f中之后，得到一个符号表达式
   gx=y(i-1)+h*fx-varvec(2); %其实这里就是方程g(x)=0左边的g(x)，就是y[n] + h*fx - y[n+1] = 0这个方程，里面y[n]是上一步知道的值，fx是含有变量y[n+1]的式子，外面的y[n+1]也视为变量
   y(i)=NewtonRoot(gx,-10,10,eps); %这就是用牛顿法来解非线性方程了

6楼: Originally posted by sudo at 2012-02-23 15:41:38:
如果它的方法要通用的话，在NewtonRoot还必须对gx进行一次符号求导（用diff函数），所以这个方法其实比较蛋疼

还不如直接用前向欧拉给估计了...

7楼: Originally posted by 微尘、梦想 at 2012-02-23 16:50:00:
嗯，我也就是想把它搞明白它到底是如何工作的，解常微分方程的显示解法，我都搞明白了，从最简单的欧拉法到四阶龙格-库塔法，不过我在网上找到一篇关于算法分析的文章，结论是隐式算法的稳定性普遍比显式算法高， ...