摘要: 先確定二分法的上下兩個端點值,這個在figure的subplot(121)中確定。
然后進行循環(huán)求解。

方程是：
 diff(s,2)+2*diff(s,1)=m*s*exp(-n*r)+h*s,
   where diff is a difference for r.
 邊界條件:
 
r=0,ds/dr=0
r=R,s=5

程序：

% 二分法+打靶法解微分方程
% 方程:
%  diff(s,2)+2*diff(s,1)=m*s*exp(-n*r)+h*s,
%    where diff is a difference for r.
% 邊界條件:
% 
% r=0,ds/dr=0
% r=R,s=5
clc;clear;close all;
m=1;
n=1;
% \copyright: zjliu
% Author's email: zjliu2001@163.com 
h=1;
R=5;
fun=inline('[s(2);m*s(1)*exp(-n*r)+h*s(1)-2*s(2)/r]',...
    'r','s','flag','m','n','h');
s_end=5;
sp=1;
rsa=0.1;
rsb=1;
d=1;
[r1,s1]=ode45(fun,[eps,R],[rsa;0],[],m,n,h);
[r2,s2]=ode45(fun,[eps,R],[rsb;0],[],m,n,h);
figure;
subplot(121);
plot(r1,s1(:,1));hold on;
plot(r2,s2(:,1),'r');
legend('for rsa','for rsb',0);
xlim([0,R]);
tg=title(['rsa=',num2str(rsa),',  rsb=',num2str(rsb)]);
subplot(122);
pa=plot(r1,s1(:,1));hold on;
pb=plot(r2,s2(:,1),'r');
pc=plot(r2,s2(:,1)*0.6,'k');
legend('for rsa','for rsb',0);
xlim([0,R]);
tg=title(['rsa=',num2str(rsa),',  rsb=',num2str(rsb)]);
while abs(s1(end,1)-s2(end,1))>1e-4; % 控制精度
    [r1,s1]=ode45(fun,[eps,R],[rsa;0],[],m,n,h);
    [r2,s2]=ode45(fun,[eps,R],[rsb;0],[],m,n,h);
    rs=[rsa+rsb]/2;
    [r,s]=ode45(fun,[eps,R],[rs;0],[],m,n,h);
    if s(end,1)>5;   % 二分法部分
        rsb=rs;
    else
        rsa=rs;
    end
    set(pc,'XData',r,'YData',s(:,1));
    set(tg,'string',['rs=',num2str(rs)]);
    pause(0.2);
end
title(['This program is over!'])