ex38ch2.m

these codes are for solving OED with matlab

function ex38ch2
x0 = 0.001;
[ys0,ysp0] = series(x0);
[x,y] = ode45(@odes,[x0 0.1],[ys0; ysp0]);
xs = linspace(x0,0.1,10);
[ys,ysp] = series(xs);
plot(x,y,xs,[ys; ysp],'o')
legend('ode45','series',2)
%=================================================
function [ys,ysp] = series(x)
a = (81/16)^(1/3);
b = (27/40)/a^2;
ys = a*x.^(4/3) + b*x.^(7/3);
ysp = (4/3)*a*x.^(1/3) + (7/3)*b*x.^(4/3);

function dydx = odes(x,y)
dydx = [ y(2); sqrt(exp(2*x)/y(1)) ];

function dfdy = Jac1(x,y)
dfdy = [  0                                        1
	     -1/2/sqrt(exp(2*x)/y(1))*exp(2*x)/y(1)^2  0 ];

function res = F2(x,y,yp)
res = [ yp(1) - y(2)
	    yp(2) - exp(x)/sqrt(y(1)) ];

function [dfdy,dfdyp] = Jac2(x,y,yp)
dfdy = [ 0                      -1
	     0.5*exp(x)*y(1)^(-3/2)  0 ];
dfdyp = [ 1   0
	      0   1 ];