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

%Find SSE for a given choice of parameters k and d.

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

%This portion sets up and solves the ODE system.

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

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);

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

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]);

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

%The following commands are used to compute and 
%add up the square errors between model and actual data.

xdata=100*y(:,2)';
actualdata = [-15.6,-14.4,-12.2,-10.4,-8.6,-6.0,-3.2,0];
z = sum((xdata-actualdata).^2);

end