求matlab高手，高分悬赏

Question

帮我看看这个程序哪里出问题了！
function [x,y]=naeuler2s(dyfun,xspan,y0,h)
%用途：改进欧拉公式解常微分方程y＇=f(x,y),y(x0)=y0
%格式：[x,y]= naeuler2(dyfun,xspan,y0,h) dyfun为函数f(x,y),xspan为求解区间[x0,xN],y0为初值y(x0),h为步长，x为返回节点，y为返回数值解
x=xspan(1):h:xspan(2);
y=zeros(length(y0),length(x));
for n=1:length(x)-1
k1= feval(dyfun,x(n),y( : ,n));
y(n+1)=y(n)+h*k1;
k2= feval(dyfun,x(n+1),y( : ,n+1));
y( : ,n+1)=y( : ,n)+h*(k1+k2)/2;
end
x=x';
y=y';
function f=dyfun(t,y)
f(1)=2*10^5*y*(log(y(1)))^2-5*10^9*y(1)*(y(1)-10^(-12));
f(2)=8*10^6-5.6*10^18*y(1)*((exp(3.63359*y(2)*10^(-9))-exp(3.63359*y(2)*10^(-9)))/2);
f=f( : );
[x,y]=naeuler2s(@dyfun,[0,1000],[1.2],0.01);
plot(x,y);

Answer 1

1、首先，需要明确一下，代码的最后两行

[x,y]=naeuler2s(@dyfun,[0,1000],[1.2],0.01);
plot(x,y);

和前面的两个函数不是同一部分，或者把函数部分保存成文件而这两行在命令行下直接运行，或者把这两行放在文件的最前面，也写成函数。

 

2、从dyfun看，所求应该是两个微分方程，那么

[x,y]=naeuler2s(@dyfun,[0,1000],[1.2],0.01);

第三个输入参数[1.2]只有一个初始条件，显然不对。

 

3、如果提供了两个初始条件，那么

f(1)=2*10^5*y*(log(y(1)))^2-5*10^9*y(1)*(y(1)-10^(-12));

等号右边有y直接参与相乘，得到结果为向量，赋值出错。你需要核实微分方程到底是什么。

 

4、目前的代码初值y0根本没有被用上。

 

5、这一行

y(n+1)=y(n)+h*k1;

疑似应为

 y(:,n+1)=y(:,n)+h*k1;

追问

我要解决的问题是这样的

初值是，y（0）=2.5*10^(-12).我不会matlab。不知道能不能帮我写一个小程序解决一下！

追答

1、z的初值？
2、你确定要求的时间范围是0-1000秒？
3、sinh有现成的函数，没必要用指数函数表示。
4、这个系统是病态的，按照你编写的欧拉法程序基本上是无法求解的，因为如果步长大了，得到的结果是错的，而步长小了，计算量又太大。