MATLAB程序出错了，怎么样修改？

Question

function y1=Q(r)
c1=2.7132e15;
M=1.989e30;
G=6.672e-11;
y1=c1^2./(r.^3)-M*G./(r.^2);
function y2=R(q)
y2=q;
function y3=S(r)
c1=2.7132e15;
y3=c1./(r.^2);
clear;clc;
format short
h=[0.05e7,0.01e7,0.005e7,0.001e7];%h为步长
for i=1:4
q=0;r(i)=0.6982e11;theta(i)=0;T=0;
while T<=50*24*3600
K1=Q(r(i)); L1=R(q);N1=S(r(i));
K2=Q(r(i)+h(i)/2*L1);L2=R(q+h(i)/2*K1);N2=S(r(i)+h(i)/2*L1);
K3=Q(r(i)+h(i)/2*L2);L3=R(q+h(i)/2*K2);N3=S(r(i)+h(i)/2*L2);
K4=Q(r(i)+h(i)*L3);L4=R(q+h(i)*K3);N4=S(r(i)+h(i)*L3);
T=T+h(i);
q=q+h(i)/6*(K1+2*K2+2*K3+K4);
r(i)=r(i)+h(i)/6*(L1+2*L2+2*L3+L4);
theta(i)=theta(i)+h(i)/6*(N1+2*N2+2*N3+N4);
end;
end
format long
disp('步长h对应取值为h=0.05e7,0.01e7,0.005e7,0.001e7时,')
disp('从远日点开始的第50天结束时水星的位置（r,theta）为')
disp(' ')
disp([r',theta'])

Answer 1

题主提供的代码没有错误，但运行时出错，其问题在于，题主把自定义函数与主程序的运行顺序搞反了。

方法1、分别将y1=Q(r)、y2=R(q)和y3=S(r)自定义函数保存为Q.m、R.m和S.m文件，再将主程序保存为main.m文件，然后在命令窗口中，执行 >>main %回车

方法2、建立main.m文件，按下列顺序书写程序

①主程序

function main( )

clear;clc;
format short 

。。。。。。

disp([r',theta'])

②Q(r)自定义函数、R(q)自定义函数和S(q)自定义函数

function y1=Q(r)
c1=2.7132e15;

。。。。。。。

function y2=R(q)
y2=q;

function y3=S(r)
c1=2.7132e15;

。。。。。。。

最后执行可得