%用打靶法求解常微分方程的边值问题，其中用到了弦割法和四阶R-K算法
function shooting_method_RK_secant
clear
clc
%close all
format long

DL=5*10^-19;  %误差判断标准，用于弦割法
XL=0.0;       %左边界点
XU=1.0;       %右边界点
h=0.01;           
N=(XU-XL)/h;
xx=XL:h:XL+N*h;
%D = 1;               %初始化误差
YL= 0.0;              %左边界点函数值
YU= 1.0;              %右边界点函数值
S0= 85;              %XL处斜率的猜测值，弦割法第一个启动点
S1= 25;              %第二个斜率猜测值，作为弦割法的第二个启动点

%-----------四阶R-K算法,将求解放在循环外部-------------------
%---------第一个试探解，求右边界的函数值及误差函数-------------------
f=@(xx,yy)[yy(2);xx^2*yy(1)+exp(-xx)]; % 定义 f1，
z10 = [YL;S0]; 
z1 = myrungekutta(f, xx, z10);
y_rk1 = z1(1,:)';   %用四阶RK算法求得的y值
z_rk1 = z1(2,:)'; 
F0=y_rk1(end)-YU; %用四阶RK算法得到的XU点的函数值与边界条件的差值
%----------第二个试探解，求右边界的函数值及误差函数--------------------
z20 = [YL;S1]; 
z2 = myrungekutta(f, xx, z20);
y_rk2 = z2(1,:)';   %用四阶RK算法求得的y值
z_rk2 = z2(2,:)'; 
F1=y_rk2(end)-YU; %用四阶RK算法得到的XU点的函数值与边界条件的差值
if (abs(F1)>abs(F0))
    y_rk=y_rk1;
else
    y_rk=y_rk2;
end
nn=0; %迭代次数
%abs(y_rk(end)-YU)
while(abs(y_rk(end)-YU)>DL)
    nn=nn+1;
    S2=S1-F1*(S1-S0)/(F1-F0); %弦割法求根
    S0=S1;
    S1=S2;
    F0=F1;
    y_rk1=y_rk;
%----------第三或更高次的试探解，求右边界的函数值及误差函数--------------------
    z20 = [YL;S1]; 
    z2 = myrungekutta(f, xx, z20);
    y_rk2 = z2(1,:)';   %用四阶RK算法求得的y值
    z_rk2 = z2(2,:)'; 
    F1=y_rk2(end)-YU %用四阶RK算法得到的XU点的函数值与边界条件的差值
    if (abs(F1)>abs(F0))
        y_rk=y_rk1;
    else
        y_rk=y_rk2;
    end
    %abs(y_rk(end)-YU)
end
hold on
plot(xx,y_rk,'r') 

function z = myrungekutta(f, xx, z0)
n = length(xx);
h = xx(2)-xx(1);
z(:,1)=z0;
for i = 1:n-1
    k1 = f(xx(i), z(:,i));
    k2=f(xx(i)+0.5*h, z(:,i) + h*k1/2);
    k3=f(xx(i)+0.5*h, z(:,i) + h*k2/2);
    k4=f(xx(i)+h, z(:,i) + h*k3);
    z(:,i+1)=z(:,i) + (k1 + 2*k2 + 2*k3 + k4)*h/6;
end

