function midsol


xx = 0:.01:1;
for J = [1 2 4 8]
  for k = 2:length(xx-1)
    x = xx(k);
    sg = exp(0.7*log(J)-2)*4*x*(1-x);
    sig(k,J) = fzero(@(sig) f1(sig,x,J),sg);
  end
end
figure(1)
plot(xx,sig(:,[1 2 4 8]))
set(gca,'Fontname','TimesNewRoman')
xlabel('x','Fontname','TimesNewRoman','Fontangle','italic')
ylabel('\sigma','Fontname','TimesNewRoman','Fontangle','italic')
h = legend('J = 1','J = 2','J = 4','J = 8');
set(h,'Fontname','TimesNewRoman')
%print -despc mid1

aa = 0.1:0.1:8;
Reg = 1e5;
for k = 1:length(aa)
  a = aa(k);
  Re(k) = fzero(@(R) f2e(R,a),Reg);
  Reg = Re(k);
end
g2e = @(a) fzero(@(R) f2e(R,a),2000)
ae = fminbnd(g2e,2,6)
Rem = g2e(ae)
aao = 0.2:0.1:16;
Rog = 5e6;
for k = 1:length(aao)
  a = aao(k);
  Ro(k) = fzero(@(R) f2o(R,a),Rog);
  Rog = Ro(k);
end
g2o = @(a) fzero(@(R) f2o(R,a),20000);
ao = fminbnd(g2o,2,6)
Rom = g2o(ao)

% numerical calculation. Vectorizing is worthwhile, as is
% FJacobian. BCJacobian seems too painful.
z = 0:0.01:1;
for bc = 1:3
  solinit = bvpinit(z,@yinit,Re(1));
  for k = 1:length(aa)
    a = aa(k)
    options = bvpset('FJacobian',@(z,y,p) fjac1(z,y,p,a),'Stats','on','Vectorized','on');
    sol = bvp5c(@(z,y,p) odefun1(z,y,p,a),@(z,y,p) bcfun1(z,y,p,a,bc),solinit,options);
    pw(bc,k) = sol.parameters;
    solinit = sol;
    sol = bvp5c(@(z,y,p) odefun1(z,y,p,a),@(z,y,p) bcfun2(z,y,p,a,bc),solinit,options);
    pt(bc,k) = sol.parameters;
    solinit = sol;
  end
end

figure(2)
Rff = (aa.^2+pi^2).^3./aa.^2;
h = plot(aa,Rff,aa,Re,aao/2,Ro/16,ae,Rem,'kx',ao/2,Rom/16,'kx',pi/sqrt(2),27/4*pi^4,'kx',aa,pw,'k.',aa,pt,'k.');
set (h,'Markersize',4)
set(gca,'Fontname','TimesNewRoman')
xlabel('a','Fontname','TimesNewRoman','Fontangle','italic')
ylabel('R','Fontname','TimesNewRoman','Fontangle','italic')
h = legend('R_{ff}','R_{rr}','R_{rf}');
set(h,'Fontname','TimesNewRoman','Fontangle','italic')
axis([0 8 0 1e4])
%print -despc mid2a
figure(4)
plot(aa,Re,aao,Ro/8,ae,Rem,'kx',ao,Rom/8,'kx')
set(gca,'Fontname','TimesNewRoman')
xlabel('a','Fontname','TimesNewRoman','Fontangle','italic')
ylabel('R','Fontname','TimesNewRoman','Fontangle','italic')
h = legend('R_e','R_o/8');
set(h,'Fontname','TimesNewRoman','Fontangle','italic')
axis([0 8 0 1e4])
%print -despc mid2b
keyboard

figure(3)
for bc = 1:3
  solinit = bvpinit(z,@yinit,Re(1));
  for k = 1:length(aa)
    a = aa(k)
    options = bvpset('FJacobian',@(z,y,p) fjac2(z,y,p,a),'Stats','on','Vectorized','on');
    sol = bvp5c(@(z,y,p) odefun2(z,y,p,a),@(z,y,p) bcfun2(z,y,p,a,bc),solinit,options);
    p3(bc,k) = sol.parameters;
    solinit = sol;
  end
end
h = plot(aa,p3,aa,Re,'--',aao/2,Ro/16,'--',aa,Rff,'--',[2.2 3.2 2.7],[763.7 2012.9 1120.9],'kx')
set(gca,'Fontname','TimesNewRoman')
xlabel('a','Fontname','TimesNewRoman','Fontangle','italic')
ylabel('R','Fontname','TimesNewRoman','Fontangle','italic')
h = legend('R_{ff}','R_{rr}','R_{rf}');
set(h,'Fontname','TimesNewRoman','Fontangle','italic')
axis([0 8 0 1e4])
save midsol


function f = f1(sig,x,J)

y = sqrt(sig+x^2);
I1p = @(z) 0.5*(besseli(0,z) + besseli(2,z));
f = 2*x^2*(x^2+y^2)*I1p(x)/besseli(0,x)*(1 - 2*x*y/(x^2+y^2)*besseli(1,x)/besseli(1,y)*I1p(y)/I1p(x)) - (x^4-y^4) - J*x*besseli(1,x)/besseli(0,x)*(1-x^2);


function f = f2e(R,a)

tau = (R/a^4).^(1/3);
q = a*sqrt(1+0.5*tau*(1+1i*sqrt(3)));
q0 = a*sqrt(tau-1);
f = imag((sqrt(3)+1i)*q.*tanh(0.5*q)) + q0.*tan(0.5*q0);


function f = f2o(R,a)

tau = (R/a^4).^(1/3);
q = a*sqrt(1+0.5*tau*(1+1i*sqrt(3)));
q0 = a*sqrt(tau-1);
f = imag((sqrt(3)+1i)*q.*coth(0.5*q)) - q0.*cot(0.5*q0);


function y = yinit(z)

s = sin(pi*z);
c = cos(pi*z);
y = [s ; pi*c ; -pi^2*s ; -pi^3*c ; pi^4*s ; pi^5*c]/(-pi^3);


function dydz = odefun1(z,y,p,a)

dydz = [y(2:6,:) ; a^6*y(1,:)-3*a^4*y(3,:)+3*a^2*y(5,:)-p*a^2*y(1,:)];


function [dfdy,dfdp] = fjac1(z,y,p,a)

dfdy = [0 1 0 0 0 0 ; 0 0 1 0 0 0 ; 0 0 0 1 0 0 ; 0 0 0 0 1 0 ; 0 0 0 0 0 1 ; a^6-p*a^2 0 -3*a^4 0 3*a^2 0];
dfdp = [0 ; 0 ; 0 ; 0 ; 0 ; -a^2*y(1)];


function res = bcfun1(ya,yb,p,a,bc)

if (bc==1) 
  % free free
  res = [ya(1) ; ya(3) ; ya(5) ; yb(1) ; yb(3) ; yb(5) ; ya(4)-1];
elseif (bc==2) 
  % rigid rigid
  res = [ya(1) ; ya(2) ; ya(5)-2*a^2*ya(3) ; yb(1) ; yb(2) ; yb(5)-2*a^2*yb(3) ; ya(4)-1];
elseif (bc==3) 
  % rigid free
  res = [ya(1) ; ya(2) ; ya(5)-2*a^2*ya(3) ; yb(1) ; yb(3) ; yb(5) ; ya(4)-1];
end


function dydz = odefun2(z,y,p,a)

dydz = [y(2:6,:) ; a^6*y(1,:)-3*a^4*y(3,:)+3*a^2*y(5,:)-3*z.^2*p*a^2.*y(1,:)];


function [dfdy,dfdp] = fjac2(z,y,p,a)

dfdy = [0 1 0 0 0 0 ; 0 0 1 0 0 0 ; 0 0 0 1 0 0 ; 0 0 0 0 1 0 ; 0 0 0 0 0 1 ; a^6-3*z^2*p*a^2 0 -3*a^4 0 3*a^2 0];
dfdp = [0 ; 0 ; 0 ; 0 ; 0 ; -3*z^2*a^2*y(1)];


function res = bcfun2(ya,yb,p,a,bc)

if (bc==1)
  % free free
  res = [ya(1) ; ya(3) ; ya(5) ; yb(1) ; yb(3) ; yb(5) ; ya(4)-1];
elseif (bc==2) 
  % rigid rigid
  res = [ya(1) ; ya(3) ; ya(4)-a^2*ya(2) ; yb(1) ; yb(3) ; yb(4)-a^2*yb(2) ; ya(4)-1];
elseif (bc==3)
  % rigid free
  res = [ya(1) ; ya(3) ; ya(4)-a^2*ya(2) ; yb(1) ; yb(3) ; yb(5) ; ya(4)-1];
end
