ex29ch3.m

these codes are for solving OED with matlab

function sol = ex29ch3
global epsilon
options = bvpset('FJacobian',@Jac,'BCJacobian',@BCJac,...
                 'Vectorized','on');
sol = bvpinit(linspace(-pi/2,pi/2,20),0.5,1);
for index = 1:10
    epsilon = 1/(10*index);
    sol = bvp4c(@ode,@bc,sol,options);
end
plot(sol.x,sol.y(1,:))
axis([-pi/2 pi/2 0 1.2]);
title(['Solution for \epsilon = ',num2str(epsilon),'.'])
fprintf('Lambda = %e.\n',sol.parameters);
%=================================================
function dydx = ode(x,y,lambda)
global epsilon
dydx = (sin(x) .^2 - lambda*sin(x) .^4 ...
         ./ y(1,:))/epsilon;

function [dFdy,dFdlambda] = Jac(x,y,lambda)
global epsilon
dFdy = lambda*sin(x)^4/y^2/epsilon;
dFdlambda = - sin(x)^4/y/epsilon;

function res = bc(ya,yb,lambda)
res = [ ya-1; yb-1 ];

function [dBCdya,dBCdyb,dBCdlambda] = BCJac(ya,yb,lambda)
dBCdya = [ 1; 0];
dBCdyb = [ 0; 1];
dBCdlambda = [ 0; 0 ];