lk2.m

希望对大家学习有所帮助

if nargin~=5
  error('wrong number of parameters in using pend2(funname,t0,step,n,y0)');
end;

y=y0; tout(1)=t0; yout(:,1)=y0;
for i=1:n
  ttemp1=tout(i);  
  ttemp2=ttemp1+step/2;
  k1=feval(funname,ttemp1,i,y);
  k2=feval(funname,ttemp2,i,y+(step/2)*k1);
  k3=feval(funname,ttemp2,i,y+(step/2)*k2);
  k4=feval(funname,ttemp1+step,i+1,y+step*k3);
  
  y=y+(step/6)*(k1+2*k2+2*k3+k4);
  yout(:,i+1)=y; 
  tout(i+1)=ttemp1+step;
end 
tout=tout';