求助四阶龙格库塔怎么出现矩阵下标错误？

comxp

function [x,y]=rk(ufunc,y0,h,a,b)
n=floor((b-a)/h); x(1)=a; y(:,1)=y0;
for ii=1:60
   x(ii+1)=x(ii)+h;
   k1=ufunc(x(ii),y(:,ii));                % 这一行出错！！
   k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);
   k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);
   k4=ufunc(x(ii)+h,y(:,ii)+h*k3);
   y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;
end
function dy=test_fun(x,y)
dy = zeros(3,1);
dy(1) = y(1) * y(3);
dy(2) = -y(1) + y(3);
dy(3) = -0.51 * y(1) * y(2);


命令窗口
>> [t,y]=rk(@test_fun,[1;1;3],0.25,0,15);
??? Subscript indices must either be real positive integers or logicals.

Error in ==> d:\MATLAB6p5\work\rk.m
On line 5  ==> k1=ufunc(x(ii),y(:,ii));

请问是什么原因？x(ii)和y(:,ii)下标应该没错啊，ii=1.

[ 本帖最后由 ChaChing 于 2009-8-8 21:11 编辑 ]

ChaChing

我用R2006a试了下, 是不会报错的! 但是plot(t,y)感觉不太对!
注意到LZ用的是v6.5, 特地进入v6.5试了下, 的确会报错!
但原因没找着, 感觉好像function handle @没成功! sorry! 这个不熟!

comxp

呵呵，难道是版本bug？一直不得其解，ls真细心，版本都看出来了~谢谢

kakalx

我在7.1下面操作；没有出现报错。但是结果不对。

ChaChing

个人水平专业有限, 不太清楚"用龙格库塔求解非线性的程序"是什含意!? 但刚刚有搜了下, 的确是有龙格库塔的程序!
个人以为用"在论坛上搜索"字眼, 并非去敷衍! 本来不就应该养成好习惯, 发问前动手搜一下, 不对吗?
再次提醒下, 没有人是有义务, 一定得帮忙自己解决问题的! 不是吗?
ps : 刚刚才发现我多挂个勋章了! 其实个人不太喜欢挂那些东东的! 原因之一就是怕许多人有这种想法!