function beertest(k,d);
  
%Michael A. Karls
%Department of Mathematical Sciences
%Ball State University
%Muncie, IN  47306 USA

%Parameters k and d are chosen as inputs for the beertest function.

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

%The function bubble is used to set up the ODE system.

function ydot=bubble(t,y,k,d);

%Enter model parameters.

g = 9.807;
P = 1.013*10^5;
R = 8.3145;
T = 281;
mCO2 = 4.41*10^(-2);
rhobeer = 1.010*10^3;
eta = 1.3*10^(-3);

%Set up the system of ODE.

ydot=[4*pi*k*((3*y(1)*R*T./(4*pi*(P-g*rhobeer*y(2)))).^2).^(1/3);
      y(3);
      g*(rhobeer*R*T./(mCO2*(P-g*rhobeer.*y(2)))-1)-d*y(3)./...
      (mCO2.*y(1))*((3*y(1)*R*T./(4*pi*(P-g*rhobeer*y(2)))).^2).^(1/3)];

end

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

%Solve the system of ODE numerically.  
%The time interval, output times, and initial data need to be specified.

tspan=[0:0.54:3.78];

x0 = -0.156;
n0 = 9.0583*10^(-10);
y0 = 0;

[t,y]=ode45(@(t,y) bubble(t,y,k,d),tspan,[n0 x0 y0]);

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

%Use the solution to the system of ODE to find bubble 
%position and radius at each time increment and display 
%bubble data in a table.

%Parameters used to compute the radius need to be redefined.

g = 9.807;
P = 1.013*10^5;
R = 8.3145;
T = 281;
rhobeer = 1.010*10^3;

disp('')

disp('time (sec), depth (cm), radius/(initial radius)')

disp([t,100*y(:,2), ...
((3*y(:,1)*R*T./(4*pi.*(P-g*rhobeer.*y(:,2)))).^(1/3))/ ...
((3*y(1,1)*R*T./(4*pi.*(P-g*rhobeer.*y(1,2)))).^(1/3))])

end