function epsilon = problem6_by_diffeq()

% approach by differential equation for generating function E(x),
% evaluated at x=1. See Maple worksheet 'GeneratingFunction.mws'
% F. Bornemann March 9, 2007

tol = 1e-14; % gives 14 correct digits in about 0.5 sec

eq = @(epsilon) ExpectedValue(epsilon) - 2;
epsilon = fzero(eq,[0.06,0.07],optimset('TolX',tol));

end

%-------------------------------------------------------------------------

function E = ExpectedValue(epsilon)

options=odeset('RelTol',100*eps,'AbsTol',eps);

global a b
a = 1/4; b = sqrt((1/4-epsilon)*(1/4+epsilon));

% integrate E(x) to x=1

[x,E]=ode113(@deq,[0,1],[1;2*(a^2+b^2)],options);

E = E(end,1);

end

%-------------------------------------------------------------------------
% differential equation for generating function E(x) from Maple:
%-------------------------------------------------------------------------

function dy = deq(x,y)

global a b;

dy = y;
if x==0
    dy(1) = y(1);
    dy(2) = 12*a^4+48*a^2*b^2+12*b^4;
else
    dy(1) = y(2);
    dy(2) = -((-24*x*b^2*a^2+12*x*b^4+12*x*a^4-2*b^2-2*a^2)*y(1) + ...
     (-96*x^2*b^2*a^2+48*x^2*a^4+48*x^2*b^4-16*x*a^2-16*x*b^2+1)*y(2))/...
     (-32*x^3*b^2*a^2+16*x^3*a^4+16*x^3*b^4-8*x^2*a^2-8*x^2*b^2+x);
end

end