打靶法求解两点边值问题

function varargout = shooting_two_point_boundary(varargin)
% ==========================================================
% 函数名：shooting_two_point_boundary.m
% 基于打靶法计算两点边值问题，仅针对二阶微分方程
% author: xianfa110.
% blog： http://blog.sina.com.cn/xianfa110
% 函数形式：
% [result,err,z0] = shooting_two_point_boundary(@fun,[y_0,y_end],[x_0,x_1],h);
% 输入：
% fun = 函数名；
% y_0 = 函数初值；
% y_end = 函数终值；
% x_0 = 自变量初值；
% x_end = 自变量终值；
% h = 积分步长；
% 输出：
% result = [x,y];
% err = 误差；
% z0 = y'初值；
% ===========================================================
% 函数fun：4y''+yy' = 2x^3 +16 ; 2<= x <=3
% 写法：
% function f = fun(y,x)
% dy = y(2);
% dz = (2*x^3+16-y(1)*y(2))/4;
% f = [dy,dz];
% ===========================================================
% 注意:y(1) = y，y(2) = y'。
% ===========================================================
F = varargin{1};
y_0 = varargin{2}(1);
y_end = varargin{2}(2);
x_0 = varargin{3}(1);
x_1 = varargin{3}(2);
ts = varargin{4};
t0 = x_0-0.5;
flg = 0;
kesi = 1e-6;
y0 = rkkt(F,[y_0,t0],x_0,x_1,ts);
n = length(y0(:,1));
if abs(y0(n,1)-y_end)<=kesi
    flg=1;
else
    t1=t0+1;
    y1=rkkt(F,[y_0,t1],x_0,x_1,ts);
    if abs(y1(n,1)-y_end)<=kesi
        flg=1;
    end
end
if flg ~= 1
    while abs(y1(n,1)-y_end) > kesi
        % ==========插值法求解非线性方程=============== %
        t2 = t1-(y1(n,1)-y_end)*(t1-t0)/(y1(n,1)-y0(n,1));
        y2 = rkkt(F,[y_0,t2],x_0,x_1,ts);
        t0=t1;
        t1=t2;
        y0=y1;
        y1=y2;
    end
end
x = x_0:ts:x_1;
out = [x',y1(:,1)];
varargout{1} = out;
varargout{2} = abs(y1(n,1)-y_end);
varargout{3} = t1;

转载：http://blog.sina.com.cn/s/blog_408540af0100b7mi.html