常微分方程组两点边值问题的数值解法
3)1(1)0(04y
y
yy
可化为微分方程组3)1(1)0(41221yyyyyy

方法一：配置法
Matlab程序：
function bvcollation
clc

solinit = bvpinit(linspace(0,1,20),[100 600]);%
sol = bvp4c(@twoode,@twobc,solinit);
x = linspace(0,1,20);
y = deval(sol,x);
y'
plot(x,y(1,:),x,y(2,:));

end
%微分方程组
function dydx = twoode(x,y)
dydx = [ y(2)
4*y(1)];
end
%边值条件
function res = twobc(ya,yb)
res = [ ya(1)-1
yb(1)-3];
end
运行结果：
1.0000 -0.4203
0.9834 -0.2117
0.9777 -0.0055
0.9828 0.2007
0.9988 0.4091
1.0259 0.6220
1.0644 0.8419
1.1147 1.0710
1.1774 1.3121
1.2531 1.5677
1.3427 1.8407
1.4472 2.1341
1.5678 2.4512
1.7057 2.7954
1.8626 3.1707
2.0401 3.5811
2.2402 4.0313
2.4652 4.5261
2.7175 5.0712
3.0000 5.6724
方法二：打靶法
程序：积分用ode45,搜索用二分法
function shoot_first_try
clc
clear all
a=0;
b=1;
a2=-1;
b2=0;
s=0;
tol=1.0*10e-6;
h1=Fun(a2)
h2=Fun(b2)
s=bisect(@Fun,a2,b2,tol)
[t1,y1]=ode45(@f,[a,b],[1,s])
plot(t1,y1(:,1),t1,y1(:,2))

function xc=bisect(f,a,b,tol)
if sign(f(a))*sign(f(b))>=0
error('f(a)*f(b)<0 not satisfied!')
end
fa=f(a);
fb=f(b);
k=0;
while (b-a)/2>tol
c=(a+b)/2;
fc=f(c);
if fc==0
break
end
if sign(f(a))*sign(f(b))<0
b=c;fb=fc;
else
a=c;fa=fc;
end
end
xc=(a+b)/2;

function z=Fun(s)
a=0;b=1;yb=3;
h=linspace(a,b,20);
y0=[1;s];
[t,y]=ode45(@f,h,y0);
z=y(end,1)-yb;

function ydot=f(t,y)
%ydot=[0;0];
ydot(1)=y(2);
ydot(2)=4*y(1);
ydot=[ydot(1);ydot(2)];

运行结果：
1.0000 -0.5000
0.9969 -0.4749
0.9940 -0.4499
0.9913 -0.4250
0.9887 -0.4001
0.9799 -0.3017
0.9736 -0.2041
0.9697 -0.1069
0.9683 -0.0100
0.9692 0.0868
0.9726 0.1839
0.9784 0.2814
0.9867 0.3797
0.9974 0.4788
1.0106 0.5792
1.0264 0.6811
1.0447 0.7846
1.0656 0.8901
1.0892 0.9978
1.1155 1.1080
1.1446 1.2210
1.1766 1.3370
1.2115 1.4564
1.2495 1.5794
1.2905 1.7064
1.3348 1.8377
1.3825 1.9735
1.4335 2.1143
1.4882 2.2603
1.5466 2.4120
1.6089 2.5698
1.6751 2.7339
1.7456 2.9049
1.8204 3.0832
1.8998 3.2692
1.9840 3.4633
2.0731 3.6661
2.1674 3.8781
2.2671 4.0998
2.3724 4.3317
2.4837 4.5745
2.5711 4.7637
2.6621 4.9596
2.7569 5.1625
2.8555 5.3726