请帮忙修改一下四阶龙格库塔法解微分方程的matlab程序

insects

请帮忙编个四阶龙格库塔法解非线性微分方程的程序给我看看，用matlab语言的，方程是
          xdot(2)=(-(1/tau1+K*tau2*(cosx(1))/tau1)*x(2)-K*(sinx(1))/tau1+Deltaomegao/tau1);
         xdot(1)=x(2);
其中xdot(1)就是x1求一阶导,xdot(2)就是x2求一阶导，Deltaomegao，tau1，tau2，K是要给定初值的常量
我把其中的cosx(1)，sinx(1)用泰勒公式展开了前三项然后用ode45编了程序，但是老师说这样不精确，也不专业，要我用龙格库塔法，不会，所以请大家帮帮忙，看看怎么帮我改改我的程序。
本人的程序如下：
function xdot=TwoPLL2(t,x,dummy,tau1,tau2,K,Deltaomegao)
         xdot=zeros(2,1);
         xdot(2)=(-(1/tau1+K*tau2*(1-x(1)^2/2+x(1)^4/24)/tau1)*x(2)-K*(x(1)-x(1)^3/6+x(1)^5/120)/tau1+Deltaomegao/tau1);
         xdot(1)=x(2);
主程序：
tau1=5;tau2=1;K=12;Deltaomegao=[1,2,3,4];
tspan=linspace(0,15,1500);
options=odeset('Abstol',[1e-6,1e-6]);
for i=1:4
     [t,x]=ode45('TwoPLL2',tspan,[-pi,Deltaomegao(i)]',options,tau1,tau2,K,Deltaomegao(i));
     figure(1);plot(t,x(:,1));       hold on; title('Phase difference as a function of Time');ylabel('θe(t)/rad');xlabel('t/s');
     figure(2);plot(t,x(:,2));       hold on; title('Frequency difference as a function of Time');ylabel('(dθe(t)/dt)/(rad?sˉ1)');xlabel('t/s');
     figure(3);plot(x(:,1),x(:,2));  hold on; title('Phase portrait');ylabel('(dθe(t)/dt)/(rad?sˉ1)');xlabel('/rad');
end:'(

insects

用matlab解微分方程组的时候，方程组中有sin,cos函数怎么办？怎么编呢？要不要先用泰勒展开先呢，ode45能不能直接解这样的方程组啊？
例：
function xdot=0000001(t,x,dummy,tau1,tau2,K,Deltaomegao)
         xdot=zeros(2,1);
         xdot(2)=(-(1/tau1+K*tau2*sin(x(1))/tau1)*x(2)-K*cos(x(1))/tau1+Deltaomegao/tau1);
         xdot(1)=x(2);
还是要展开先呢？将下程序保存为TwoPLL2.m
function xdot=TwoPLL2(t,x,dummy,tau1,tau2,K,Deltaomegao)
         xdot=zeros(2,1);
         xdot(2)=(-(1/tau1+K*tau2*(1-x(1)^2/2+x(1)^4/24)/tau1)*x(2)-K*(x(1)-x(1)^3/6+x(1)^5/120)/tau1+Deltaomegao/tau1);
         xdot(1)=x(2);
不过自己展开的话精确度就小了很多，因为自己展开不有写出很多项，只能展开前几项啊，能解决这个问题吗？
主程序：
tau1=5;tau2=1;K=12;Deltaomegao=[1,2,3,4,5,6,7,8,9,10];
tspan=linspace(0,15,1500);
options=odeset('Abstol',[1e-6,1e-6]);
for i=1:10
     [t,x]=ode15s('TwoPLL2',tspan,[-pi,Deltaomegao(i)]',options,tau1,tau2,K,Deltaomegao(i));
     figure(1);plot(t,x(:,1));       hold on; title('Phase difference as a function of Time');ylabel('θe(t)/rad');xlabel('t/s');
     figure(2);plot(t,x(:,2));       hold on; title('Frequency difference as a function of Time');ylabel('(dθe(t)/dt)/(rad?sˉ1)');xlabel('t/s');
     figure(3);plot(x(:,1),x(:,2));  hold on; title('Phase portrait');ylabel('(dθe(t)/dt)/(rad?sˉ1)');xlabel('/rad');
end
可以运行一下，没有错误的，我有运行过了的。:handshake

insects

将下程序保存为TwoPLL2.m
function xdot=TwoPLL2(t,x,dummy,tau1,tau2,K,Deltaomegao)
         xdot=zeros(2,1);
         xdot(2)=(-(1/tau1+K*tau2*(1-x(1)^2/2+x(1)^4/24)/tau1)*x(2)-K*(x(1)-x(1)^3/6+x(1)^5/120)/tau1+Deltaomegao/tau1);
         xdot(1)=x(2);
程序ode4.m  四阶龙格库塔法
function Y = ode4(odefun,tspan,y0,varargin)
if ~isnumeric(tspan)
error('TSPAN should be a vector of integration steps.');
end
if ~isnumeric(y0)
error('Y0 should be a vector of initial conditions.');
end
h = diff(tspan);
if any(sign(h(1))*h <= 0)
error('Entries of TSPAN are not in order.')
end
try
f0 = feval(odefun,tspan(1),y0,varargin{:});
catch
msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];
error(msg);
end
y0 = y0(:); % Make a column vector.
if ~isequal(size(y0),size(f0))
error('Inconsistent sizes of Y0 and f(t0,y0).');
end
neq = length(y0);
N = length(tspan);
Y = zeros(neq,N);
F = zeros(neq,4);
Y(:,1) = y0;
for i = 2:N
ti = tspan(i-1);
hi = h(i-1);
yi = Y(:,i-1);
F(:,1) = feval(odefun,ti,yi,varargin{:});
F(:,2) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,1),varargin{:});
F(:,3) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,2),varargin{:});
F(:,4) = feval(odefun,tspan(i),yi+hi*F(:,3),varargin{:});
Y(:,i) = yi + (hi/6)*(F(:,1) + 2*F(:,2) + 2*F(:,3) + F(:,4));
end
Y = Y.';
主程序：
tau1=5;tau2=1;K=12;Deltaomegao=[1,2,3,4,5,6,7,8,9,10];
tspan=linspace(0,15,1500);
options=odeset('Abstol',[1e-6,1e-6]);
for i=1:10
     [t,x]=ode4('TwoPLL2',tspan,[-pi,Deltaomegao(i)]',options,tau1,tau2,K,Deltaomegao(i));
     figure(1);plot(t,x(:,1));       hold on; title('Phase difference as a function of Time');ylabel('θe(t)/rad');xlabel('t/s');
     figure(2);plot(t,x(:,2));       hold on; title('Frequency difference as a function of Time');ylabel('(dθe(t)/dt)/(rad?sˉ1)');xlabel('t/s');
     figure(3);plot(x(:,1),x(:,2));  hold on; title('Phase portrait');ylabel('(dθe(t)/dt)/(rad?sˉ1)');xlabel('/rad');
end
本人找到一个四阶龙格库塔法的程序，但是在主程序中调用时出现如下错误，:@L 请教一下怎么修改？:@L
??? Error using ==> ode4
Too many output arguments.
Error in ==> TwoPLLshzhjiyituduxi at 5
     [t,x]=ode4('TwoPLL2',tspan,[-pi,Deltaomegao(i)]',options,tau1,tau2,K,Deltaomegao(i));