lorenz混沌系统 MATLAB仿真

function [t,y] = odeRK4sysv(diffeq,tn,h,y0,varargin) t = (0:h:tn)'; nt = length(t); neq = length(y0); y = zeros(nt,neq); y(1,:) = y0(:)'; flag = []; h2 = h/2; h3 = h/3; h6 = h/6; k1 = zeros(neq,1); k2 = k1; k3 = k1; k4 = k1; ytemp = k1; for j=2:nt told = t(j-1); yold = y(j-1,:)'; k1 = feval(diffeq,told,yold,flag,varargin{:}); ytemp = yold + h2*k1; k2 = feval(diffeq,told+h2,ytemp,flag,varargin{:}); ytemp = yold + h2*k2; k3 = feval(diffeq,told+h2,ytemp,flag,varargin{:}); ytemp = yold + h*k3; k4 = feval(diffeq,told+h,ytemp,flag,varargin{:}); y(j,:) = ( yold + h6*(k1+k4) + h3*(k2+k3) )'; end