ex21ch3.m

these codes are for solving OED with matlab

function ex21ch3
solinit = bvpinit(linspace(0,1,10),@guess,0.5);
sol = bvp4c(@odes,@bcs,solinit);
t0 = sol.parameters;
fprintf('t0 is about %7.4f.\n',t0)
fprintf('The reduced drag coefficient is about %7.4f.\n',sol.y(3,1))
t = t0 + (1 - t0)*sol.x;
t = [0 t];
y = [0 sol.y(1,:)];
plot(t,y)
figure
theta = linspace(0,2*pi,20);
X = t' * cos(theta);
Y = t' * sin(theta);
Z = y' * ones(1,20);
surf(X,Y,Z)
view(-45,-40)
grid off
%=============================================
function v = guess(x)
v = [ x; 1; 0.5^2 ];

function dydx = odes(x,y,t0)
t = t0 + (1 - t0)*x;
dydx = [ y(2)
         y(2)*(1 + y(2)^2)/(t*(3*y(2)^2 - 1))
         2*t - 2*t/(1 + y(2)^2) ];
dydx = (1 - t0)*dydx;

function res = bcs(ya,yb,t0)
res = [ ya(1); ya(2)-1; yb(1)-1; yb(3)-1 ];