Iteration to convergence for a steam engine system

Hello ladies and gentlemen of the matlab community. I am currently working on an assignment that requires all the forces acting within a piston steam engine to be calculated and graphed using matlab. I am attempting to iterate a system for thetadot until the resulting Kinetic Energy is equal to what it should. Please find below my buggy code. I would be terribly appreciative to anyone who can fix my code asap :)

R = 0.699;  
theta = 3;  
thetaDot = 500*(2*pi/60);  
Lrod = 2*R + 0.06;
RodMass = 3.4257;
PistonMass = 0.6176;
IRod = 0.6195;
IFly = 67.25;
Beta = asin((R*sin(theta))/Lrod);
x = (3*R + 60) - sqrt((Lrod^2)+(R^2)+Lrod*R*cos(theta+Beta));
OA = 0.010800000000000;
OB = 0.0258;
OC = 1.398;

%define betadot
betaDot = ((R*cos(theta)*thetaDot)/(Lrod*cos(Beta)))*thetaDot;

%define xdot
xDot = ((Lrod^2 + R^2 + Lrod*R*cos(theta+Beta))^(-1/2))*(Lrod+2*R*sin(theta+Beta)*(thetaDot+betaDot))*thetaDot;

%KE of piston
KEPist = (1/2)*PistonMass*xDot^2;

%Vrod?

%KE of rod
KERod = (1/2*RodMass*(xDot*thetaDot)^2)+(1/2*IRod*betaDot^2);

%KE of flywheel
KEFly = (1/2)*IFly*thetaDot^2;

KEi = 0.5*IFly*(500*(2*pi*(R+0.03)/60))^2;

if x < OB && theta < pi
    deltaKE1 = WorkinAB - Workamb;
elseif x >= OB && x <= OC
    deltaKE1 = WorkinBC - Workamb;
elseif theta > pi;
    deltaKE1 = 'something';

end

deltaKE = 100;

KE1 = KEi + deltaKE;

deltaThetaDot = 50;

lim = 10;

KE2 = 0;

count = 0;
sign1 = 1;

while abs(deltaKE) > lim;

      if KE1 > KE2

          thetaDot = thetaDot + deltaThetaDot;

      else

          thetaDot = thetaDot - deltaThetaDot;
      end

      KE2 = KEPist + KEFly + KERod;

      deltaKE = KE1 - KE2;

      if deltaKE > 0
          sign2 = 1;
      else
          sign2 = 0;
      end

      if sign2 ~= sign1
          deltaThetaDot = (deltaThetaDot)/2;
      end

      sign1 = sign2;

      %define betadot
  betaDot = ((R*cos(theta)*thetaDot)/(Lrod*cos(Beta)))*thetaDot;

%define xdot
xDot = ((Lrod^2 + R^2 + Lrod*R*cos(theta+Beta))^(-1/2))*(Lrod+2*R*sin(theta+Beta)*(thetaDot+betaDot))*thetaDot;

%KE of piston
KEPist = (1/2)*PistonMass*xDot^2;

%Vrod?

%KE of rod
KERod = (1/2*RodMass*(xDot*thetaDot)^2)+(1/2*IRod*betaDot^2);

%KE of flywheel
KEFly = (1/2)*IFly*thetaDot^2;

%dKE2 (difference between work in and work out
KE2 = KEPist + KERod + KEFly;

count = count+1;
disp(count);
disp(deltaKE);
disp(thetaDot);

if count == 100;
    return
end
end

[Walter Roberson]

What error message are you encountering, or what difference do you see between what you are getting and what you expect to get?