diff --git a/SIRC.m b/SIRC.m
index 6adb0ba8d3b51195941e68eb71620821c3c7f7d1..0ffd14dd9ae3dd9373e3e2ec665f873fffb7abf7 100644
--- a/SIRC.m
+++ b/SIRC.m
@@ -11,7 +11,7 @@ C0 = 0;
 S0 = N - I0 - R0 - C0;
 % Initial parameteres
 nu = 14;
-mu = 0.015;
+mu = 0.013;
 epsilon = 1.2;
 beta = 0.0002;
 gamma = 0.05;
@@ -209,5 +209,14 @@ ylabel('Number of individuals');
 legend('S', 'I', 'R', 'C');
 
 %% Gradient descent method
-% For dCdt = gamma*q*I - (Gamma-mu)*C
-theta_C = gd(X_C2, Y_C2, [gamma_lsq2*q_lsq2, Gamma_lsq2-mu_lsq2]', 0.01, 5);
+
+% dSdt + dIdt = nu - mu*S - (gamma+mu)*I
+theta_SI = abs(gd(X_SI2, Y_SI2, [nu_lsq2, mu_lsq2, gamma_lsq2+mu_lsq2]', 0.01, 1));
+% dIdt = beta*I*S + epsilon*beta*C*S - (gamma+mu)*I
+theta_I = abs(gd(X_I2, Y_I2, [beta_lsq2, epsilon_lsq2*beta_lsq2, gamma_lsq2+mu_lsq2]', 0.01, 1));
+% dCdt = gamma*q*I - (Gamma-mu)*C
+theta_C = abs(gd(X_C2, Y_C2, [gamma_lsq2*q_lsq2, Gamma_lsq2-mu_lsq2]', 0.01, 1));
+
+[nu_gd, mu_gd, beta_gd, epsilon_gd, gamma_gd, q_gd, Gamma_gd] = deal(theta_SI(1), ...
+    theta_SI(2), theta_I(1), theta_I(2)/theta_I(1), theta_SI(3)-theta_SI(2), ...
+    theta_C(1)/(theta_SI(3)-theta_SI(2)), theta_C(2)+theta_SI(2));
\ No newline at end of file
diff --git a/gd.m b/gd.m
index c5b1088be900778da34fe6c6d3be1c8eef6b4f0b..205516f22272f8134b2fdf06a905b54a3b0cc5f6 100644
--- a/gd.m
+++ b/gd.m
@@ -2,6 +2,6 @@ function theta = gd(X, y, theta, alpha, num_iters)
 m = length(y);
 for iter = 1:num_iters
      for i=1:size(X,2)
-         theta(i) = theta(i) - alpha/m * sum((X*theta - y).*X(:,i));
+         theta(i) = theta(i) - alpha/m * sum((X*theta-y)).*X(iter, i);
      end
  end
\ No newline at end of file