function [z,y] = sol2(a)

format long
options = odeset('Events',@events,'RelTol',1e-12,'AbsTol',1e-6);

figure(1)
[z,y] = ode45(@plume1,[0 3], [eps eps 1],options);
z(end)
q0 = 0.6*(9/20)^(1/3);
m0 = (9/20)^(2/3);
b0 = 1.5*q0;
y0 = [q0*eps^(5/3) m0*eps^(4/3) 1+b0*eps^(8/3)]
[z2,y2] = ode45(@plume1,[eps 3], y0,options);
z2(end)
plot(y,z,y2,z2)
legend('Q','M','B','Location','SouthWest')
axis([min(y(:,3)) max(y(:,1)) 0 z(end)])

figure(2)
[z,y] = ode45(@plume2,[0 3], [0 0 1],options);
z(end)
Q = y(:,1).*z.^(2/3);
M = y(:,2).*z.^(1/3);
B = y(:,3);
plot(Q,z,M,z,B,z)
legend('Q','M','B','Location','SouthWest')
axis([min(B) max(Q) 0 z(end)])

figure(3)
[z,y] = ode45(@(z,y) plume3(z,y,a),[0 100], [0 0 0],options);
z(end)
Q = y(:,1).*z.^(2/3);
M = y(:,2).*z.^(1/3);
B = 1 + y(:,3).*z.^(a+5/3);
plot(Q,z,M,z,B,z)
legend('Q','M','B','Location','SouthWest')
%axis([min(B) max(Q) 0 z(end)])


function yp = plume1(z,y)

Q = y(1); M = y(2); B = y(3);
yp = [abs(M)^0.5 ; B*Q/M ; -Q];


function yp = plume2(z,y)

if (z==0)
  q0 = 0.6*(9/20)^(1/3);
  m0 = (9/20)^(2/3);
  yp = [q0 ; m0 ; 0];
else
  q = y(1); m = y(2); B = y(3);
  yp = [z^(-0.5)*abs(m)^0.5 - 2/3/z*q ; B*q/m - 1/3/z*m ; -z^(2/3)*q];
end


function yp = plume3(z,y,a)

if (z==0)
  q0 = 0.6*(9/20)^(1/3);
  m0 = (9/20)^(2/3);
  b0 = -q0/(a-2/3);
  yp = [q0 ; m0 ; b0];
else
  q = y(1); m = y(2); b = y(3); B = 1 + b*z^(a+5/3);
  yp = [z^(-0.5)*abs(m)^0.5 - 2/3/z*q ; B*q/m - 1/3/z*m ; -1/z*q - (a+5/3)/z*b];
end


function [value,isterminal,direction] = events(t,y)

value = y(2);     % detect height = 0
isterminal = 1;   % stop the integration
direction = -1;   % negative direction


