www.gusucode.com > funfun工具箱matlab源码程序 > funfun/private/ntrp4h.m

    function [yint,ypint] = ntrp4h(tint,t,y,tnew,ynew,ymid,yp,ypnew)
%NTRP4H  Interpolation helper function for BVP5C.
%   YINT = NTRP4H(TINT,T,Y,TNEW,YNEW,YMID,YP,YPNEW) evaluates the quartic
%   interpolant interpolant at time TINT. TINT may be a scalar or a row
%   vector.  The quartic interpolates y,yp at t; ynew,ypnew at tnew; and
%   ymid at (t+tnew)/2.
%   [YINT,YPINT] = NTRP4H(TINT,T,Y,TNEW,YNEW,YMID,YP,YPNEW) returns also 
%   the derivative of the interpolating polynomial. 
%   
%   See also BVP5C, DEVAL.

%   Jacek Kierzenka and Lawrence F. Shampine
%   Copyright 1984-2007 The MathWorks, Inc.

% Convert to the scaled variable s with x = t + sh.  Must then convert 
% the derivatives: d/ds = d/dx * dx/ds = h*d/dx.
h = tnew - t;
s = (tint - t)/h;
s2 = s .* s;
s3 = s .* s2;
s4 = s .* s3;
y0p = h*yp;
y1p = h*ypnew;
del1 = ymid - y;
del2 = ynew - ymid;
del3 = y1p - y0p;
A2 = ((  11*del1 -  5*del2) +   del3) - 3*y0p;
A3 = ((- 18*del1 + 14*del2) - 3*del3) + 2*y0p;
A4 =  (   8*del1 -  8*del2) + 2*del3;
yint = y(:,ones(size(tint))) + (y0p*s + A2*s2 + A3*s3 + A4*s4);
if nargout > 1
  ypint = y0p(:,ones(size(tint))) + (2*A2*s + 3*A3*s2 + 4*A4*s3); 
  % Convert from d/ds to d/dx:  
  ypint = ypint ./ h;
end