diff --git a/CavitationPID.m b/CavitationPID.m
new file mode 100644
index 0000000..5016d25
--- /dev/null
+++ b/CavitationPID.m
@@ -0,0 +1,90 @@
+global error radius_d dt time_array tex MaxRadExp;
+
+error = 0;
+        
+% Reading experimental data
+% data2 = bubbleRadiuscell{i,j}';
+% data1 = bubbleTimeVectorcell{i,j}';
+
+
+% Find the maximum radius and the corresponding time of maximum occurence
+[MaxRadExp,I_MaxRadExp] = max(data2);
+
+
+tex = data1(I_MaxRadExp:end)-data1(I_MaxRadExp);
+rex = data2(I_MaxRadExp:end);
+
+radius_d = rex;
+dt = 0.00001;
+time_array = 0:dt:dt*(numel(tex)-1);
+
+%% CONSTANTS %%
+R_0 = MaxRadExp;
+% R_0 = 1.03813e-03;
+mu = 1e-3;
+n = 3;
+D = 1.76e-09;
+K_B = 6.737e09;
+alpha_M = .8544;
+R_v = 461.52;
+Na = 6.02e23;
+Rg = 8.314;
+T_infty = 293.15; 
+P_infty = 1.01e05; 
+Pinf = 1.01e05; 
+T_0 = T_infty;
+TB_0 = T_infty;
+Vol_0 = 4/3*pi*R_0^3;
+P_0 = 50000; %%%%%%%%%%%%%%%%%%%
+P_water = 0.83775*55000; %%%%%%%%%%%%%%%%%%%
+n_t_0 = P_0*Vol_0*Na/(Rg*T_0);
+n_water_0 = P_water / P_0 * n_t_0;
+n_air_0 = n_t_0 - n_water_0;
+
+% Van der Waals Constants
+a_air_VdW = 0.1402;
+a_air_water_VdW = 0.278594;
+a_water_VdW = 0.5536;
+b_air_VdW = 3.753e-05;
+b_air_water_VdW = 3.38882e-05;
+b_water_VdW = 3.049e-05;
+
+% Initial Conditions
+R_0 = R_0;
+Rdot_0 = 0;
+mdot_0 = 0;
+n_t_0 = P_0*Vol_0*Na/(Rg*T_0);
+n_water_0 = P_water / P_0 * n_t_0;
+n_air_0 = n_t_0 - n_water_0;
+a_0_VdW = a_air_VdW*(n_air_0/n_t_0)^2 + 2*a_air_water_VdW*( n_air_0* n_water_0/n_t_0^2) + a_water_VdW*(n_water_0/n_t_0)^2;
+b_0_VdW = b_air_VdW*(n_air_0/n_t_0)^2 + 2*b_air_water_VdW*( n_air_0* n_water_0/n_t_0^2) + b_water_VdW*(n_water_0/n_t_0)^2;
+T_0 = T_infty;
+P_0 = 50000; 
+Vol_0 = 4/3*pi*R_0^3;
+TB_0 = T_infty;
+
+% Initialization
+x0 = [n_air_0; n_water_0; n_t_0; a_0_VdW; b_0_VdW; T_0; Vol_0; TB_0; P_0; mdot_0; R_0; Rdot_0];
+
+
+%% Main Section: Solving the Equations
+tic
+options = odeset('AbsTol',1e-3,'Stats','on');%'Reltol',1e-9,
+[t_ode, x] = ode45(@ODESolver_PID, [time_array(1),time_array(end)], x0, options);
+
+% Error calculation
+radius_desired_ode = interp1(tex,radius_d,t_ode);
+e = x(:,11) - radius_desired_ode; % Error radius
+toc
+
+figure(205)
+plot(t_ode,x(:,11),'-r','LineWidth',2.75)
+hold on 
+h=plot(tex,rex,'o', 'MarkerSize',8,'Color',[0, 0.4470, 0.7410]),set(h, 'markerfacecolor', get(h, 'color'));
+title('Bubble Radius','fontsize',18,'FontName','Garamond','FontWeight','bold');
+xlabel('time (sec)', 'fontsize',28,'FontName','Garamond','FontWeight','bold');
+ylabel('Radius (m)', 'fontsize',28,'FontName','Garamond','FontWeight','bold');
+legend('Predicted Radius (MODEL)', 'Desired Radius (EXP)')
+grid on
+set(gcf,'Color',[1 1 1]);
+set(gca,'fontsize',28,'FontName','Garamond','FontWeight','bold','Color','w');
\ No newline at end of file
diff --git a/ODESolver_PID.m b/ODESolver_PID.m
new file mode 100644
index 0000000..aed7fb5
--- /dev/null
+++ b/ODESolver_PID.m
@@ -0,0 +1,146 @@
+function dx = ODESolver_PID(t, x)
+persistent radius_dPrev time_prev 
+
+if isempty(radius_dPrev)
+    radius_dPrev = 0;
+    time_prev = 0;
+end
+
+global error radius_d tex MaxRadExp; 
+
+dx = zeros(12,1);
+
+%Parameters:
+sigma = 0.072;                  % Surface tension, kg/s^2
+mu = 1e-3;                      % Water kinematic viscosity, m^2/s
+rho_L = 1e3;                    % Water density, kg/m^3
+k = 1.4;                        % Specific heat ratio
+c = 1500;
+R_0 = MaxRadExp;
+n = 3;
+D = 1.76e-09;
+K_B = 6.737e09;
+alpha_M = .8544;
+R_v = 461.52;
+Na = 6.02e23;
+Rg = 8.314;
+Rdot_0 = 0;
+mdot_0 = 0;
+T_infty = 293.15; 
+P_infty = 1.01e05; 
+Pinf = 1.01e05; 
+T_0 = T_infty;
+TB_0 = T_infty;
+Vol_0 = 4/3*pi*R_0^3;
+P_0 = 50000; 
+P_water = 0.83775*55000; 
+n_t_0 = P_0*Vol_0*Na/(Rg*T_0);
+n_water_0 = P_water / P_0 * n_t_0;
+n_air_0 = n_t_0 - n_water_0;
+% Van der Waals Constants
+a_air_VdW = 0.1402;
+a_air_water_VdW = 0.278594;
+a_water_VdW = 0.5536;
+b_air_VdW = 3.753e-05;
+b_air_water_VdW = 3.38882e-05;
+b_water_VdW = 3.049e-05;
+a_0_VdW = a_air_VdW*(n_air_0/n_t_0)^2 + 2*a_air_water_VdW*( n_air_0* n_water_0/n_t_0^2) + a_water_VdW*(n_water_0/n_t_0)^2;
+b_0_VdW = b_air_VdW*(n_air_0/n_t_0)^2 + 2*b_air_water_VdW*( n_air_0* n_water_0/n_t_0^2) + b_water_VdW*(n_water_0/n_t_0)^2;
+
+if (x(6)<258.15)
+    P_v = 180.547;
+elseif (x(6)>258.15 && x(6)<647.3)
+    P_v = 133.322*exp(18.3036-3816.44/(x(6)-46.13));
+else
+    P_v = 2.075e07;
+end
+P_v_infty = 133.322*exp(18.3036-3816.44/(T_infty-46.13));
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%% PID IMPLEMENTATION %%%%%%%%%%%%%%%%%%%%%%%%%% 
+% radius_d first derivative
+radius_desired_current = interp1(tex,radius_d,t);
+if t==time_prev
+    radius_dDot = 0;
+else
+    radius_dDot  = ( radius_desired_current - radius_dPrev ) / (t-time_prev);
+end
+radius_dPrev = radius_desired_current;
+
+% PID tuning
+Kp = 6500000;
+Kd = 0;%25;%20;
+Ki = 0;%250000;
+
+ind = find(tex<t & tex>time_prev);
+if ~isempty(ind) 
+
+    % u: feedback signal
+    u = Kp*(radius_desired_current - x(11)) + Kd*( radius_dDot - x(12)) + Ki*error;
+    error = error + (radius_desired_current - x(1))*(t-time_prev);
+else
+    u = (Kp*(radius_desired_current - x(11)) + Kd*( radius_dDot - x(12)) + Ki*error)/2.5;
+    error = error + (radius_desired_current - x(1))*(t-time_prev);
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%% Equations to solve %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Y = x(5)/(Na*x(7)/x(3)) + 0.625*(x(5)/(Na*x(7)/x(3)))^2 + 0.2869*(x(5)/(Na*x(7)/x(3)))^3+0.115*(x(5)/(Na*x(7)/x(3)))^4;
+kappa0_water = 9.98e-05 * x(8) - 0.0119;
+kappa0_air = 5.39e-05 * x(8) + 0.0108;
+kappa0_mix = 0.5*(x(2)/x(3)*kappa0_water + x(1)/x(3)*kappa0_air + 1/(x(2)/(x(3)*kappa0_water)+x(1)/(x(3)*kappa0_air)));
+kappa = x(5)/(Na*x(7)/x(3))*(1/Y+1.2+0.755*Y)*kappa0_mix;
+
+if (0.1 * P_infty - x(1) / x(3) * x(9) <0)
+	Z = x(8);
+else 
+	Z = T_0;
+end
+% keyboard;
+G = -1/(2*1.138e-23) * sqrt(pi/1.38e-23) * 1.173 * kappa * 4e-19 / n;
+lambda = x(7) / (sqrt(2)* 0.4e-18 * x(3));
+
+E_water_infty = 1.04334e-20 + 6.27210e-23 * (T_infty-250);   
+E_water_T_B = 1.04334e-20 + 6.27210e-23 * (x(8)-250);   
+E_air_Z = 8.64814e-21 + 4.40808e-23 * (Z-250);   
+E_air_T = 8.64814e-21 + 4.40808e-23 * (x(6)-250);   
+E_water_T = 1.04334e-20 + 6.27210e-23 * (x(6)-250);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ODEs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+dx(1) = 4*pi*x(11)*D*1000*rho_L*Na/(18*K_B)*(0.1*P_infty-x(1)/x(3)*x(9));
+
+GammaFunc = x(10)*x(3)/(x(2)*x(9))*sqrt(R_v*x(8)/2);
+funGamma = @(ttt) exp(-ttt.^2);
+IntTerm = integral(@(ttt)funGamma(ttt),0, GammaFunc);
+Gamma = exp(-GammaFunc^2)-sqrt(pi)* GammaFunc*(1-2/sqrt(pi)*IntTerm);
+
+dx(2) = 4*pi*x(11)^2*1000*Na/18*alpha_M/sqrt(2*pi*R_v)*(P_v_infty/sqrt(T_infty)- Gamma*x(2)*x(9)/(x(3)*sqrt(x(8))));
+dx(3) = dx(1)+ dx(2);
+
+M = 0.001/(Na*x(3))*(18*x(2)+28.97*x(1));
+dM = 0.001/(Na*x(3)^2)*(18*(x(3)*dx(2)-x(2)*dx(3))+28.97*(x(3)*dx(1) -x(1)*dx(3)));
+X = Na*x(7)/x(3);
+dX = 4*pi/3*Na/x(3)^2*(3*x(11)^2*x(12)*x(3)-x(11)^3*dx(3));
+
+dx(4) = 2*a_air_VdW*x(1)/x(3)^3*(x(3)*dx(1)-x(1)*dx(3)) + 2*a_air_water_VdW*x(2)/x(3)^3*(x(3)*dx(1)-x(1)*dx(3)) + 2*a_air_water_VdW*x(1)/x(3)^3*(x(3)*dx(2)-x(2)*dx(3))+ 2*a_water_VdW*x(2)/x(3)^3*(x(3)*dx(2)-x(2)*dx(3));
+dx(5) = 2/x(3)^3*((b_air_VdW*x(1)+b_air_water_VdW*x(2))*(x(3)*dx(1)-x(1)*dx(3))+( b_water_VdW*x(2)+b_air_water_VdW*x(1))*(x(3)*dx(2)-x(2)*dx(3)));
+dx(6) = 1/(4.40808e-23*x(1)+6.2721e-23*x(2))*(-x(9)*4*pi*x(11)^2*x(12)+kappa*4*pi*x(11)^2*(x(8)-x(6))/(n*lambda)+4*pi*x(11)^2*1000/18*Na*alpha_M/(sqrt(2*pi*R_v))*(P_v_infty/sqrt(T_infty)*E_water_infty- Gamma*x(2)*x(9)/(x(3)*sqrt(x(8)))*E_water_T_B)+4*pi*x(11)*D*1000*rho_L*Na/(18*K_B)*(0.1*P_infty-x(1)/x(3)*x(9))*E_air_Z-x(3)^2*x(4)/Na^2/x(7)^2*4*pi*x(11)^2*x(12)+1/(Na^2*x(7))*(2*x(3)*dx(3)*x(4)+x(3)^2*dx(4))-dx(1)*E_air_T-dx(2)*E_water_T);
+dx(7) = 4*pi*x(11)^2*x(12);
+dx(8) = G*(-dx(6)* sqrt(M/x(8))+(x(8)-x(6))/(2*x(8))* sqrt(x(8)/M)*dM)/(1-G*sqrt(M/x(8))+M/(2*x(8)^2)*(x(8)-x(6))*sqrt(x(8)/M));
+dx(9) = (Rg*dx(6)-(x(9)+x(4)/X^2)*(dX-dx(5)))/(X-x(5))+2*x(4)/X^3*dX-dx(4)/X^2;
+dx(10) = -alpha_M/sqrt(2*pi*R_v)*Gamma*(x(9)/sqrt(x(8))/x(3)^2*(x(3)*dx(2)-x(2)*dx(3))+x(2)/x(3)*(dx(9)/sqrt(x(8))-0.5*x(9)/x(8)^1.5*dx(8)));
+
+Rho_g = 0.001/(Na*x(7))*(18*x(2)+28.97*x(1));
+dRho_g = 0.001/(Na*x(7)^2)*(x(7)*(18*dx(2)+28.97*dx(1))-dx(7)*(18*x(2)+28.97*x(1)));
+
+
+dx(11) = x(12);
+
+% dx(12) = ((1+x(12)/c)/rho_L*(x(9)-2*sigma/x(11)-4*mu/x(11)*(x(12)-x(10)/rho_L)-x(10)^2*(1/rho_L-1/Rho_g)-Pinf)+dx(10)*x(11)/rho_L*(1-x(12)/c+x(10)/(rho_L*c))+ x(10)/rho_L*(x(12)+ x(10)/(2*rho_L)+x(12)* x(10)/(2*rho_L*c))-1.5*x(12)^2*(1-x(12)/(3*c)+2* x(10)/(3*rho_L*c))+x(11)/(rho_L*c)*(dx(9)+2*sigma*x(12)/x(11)^2+4*mu*(x(12)^2/x(11)^2+(x(11)*dx(10)-x(12)* x(10))/(rho_L*x(11)^2))-2*x(10)*dx(10)*(1/rho_L-1/Rho_g)-x(10)^2/Rho_g^2*dRho_g))/(x(11)*(1-x(12)/c+x(10)/(c*rho_L))+4*mu/(c*rho_L));
+%%%%%%%%%%%%%%%%%%%%%%%%%%% ADDING FEEDBACK SIGNAL %%%%%%%%%%%%%%%%%%%%%%%%%%%
+dx(12) = (u+(1+x(12)/c)/rho_L*(x(9)-2*sigma/x(11)-4*mu/x(11)*(x(12)-x(10)/rho_L)-x(10)^2*(1/rho_L-1/Rho_g)-Pinf)+dx(10)*x(11)/rho_L*(1-x(12)/c+x(10)/(rho_L*c))+ x(10)/rho_L*(x(12)+ x(10)/(2*rho_L)+x(12)* x(10)/(2*rho_L*c))-1.5*x(12)^2*(1-x(12)/(3*c)+2* x(10)/(3*rho_L*c))+x(11)/(rho_L*c)*(dx(9)+2*sigma*x(12)/x(11)^2+4*mu*(x(12)^2/x(11)^2+(x(11)*dx(10)-x(12)* x(10))/(rho_L*x(11)^2))-2*x(10)*dx(10)*(1/rho_L-1/Rho_g)-x(10)^2/Rho_g^2*dRho_g))/(x(11)*(1-x(12)/c+x(10)/(c*rho_L))+4*mu/(c*rho_L));
+
+time_prev = t;
+end
\ No newline at end of file