syms x1 x2 x3 x4
A = ((x1*x4) + (x2*x4) - (x1*x3));
I = ((x1*(x4)^3 + x2*(x4)^3 - x1*(x3)^3)/(12));
L = 40;
P_load = 10;
rho = 0.29;
E = 30*10^6;
w = P_load + rho*A;
W = rho*L*A;
delta = (w*(L^4))/(8*E*I);
sigma_m = (w*(L^2)*x4)/(4*I);
freq = (3.52/(2*pi))*(sqrt((E*I*386.4)/(rho*A*L^4)));
f = W-10;
g2 = (delta)-.1;
g3 = (freq)-80;
g = [g2;g3];
grad_f = gradient(f,[x1 x2 x3 x4]);
grad_g = [gradient(g2,[x1 x2 x3 x4]) gradient(g3,[x1 x2 x3 x4])];
g_x0 = vpa(subs(g,[x1 x2 x3 x4],[0.2 0.2 2.8 3.6]),4);
grad_g_x0 = vpa(subs(grad_g,[x1 x2 x3 x4],[0.2 0.2 2.8 3.6]),4);
Q = vpa((grad_g_x0*((grad_g_x0'*grad_g_x0)^-1)),4);
s2 = vpa(-Q*g_x0,4);
x_new = vpa([0.2+s2(1); 0.2+s2(2); 2.8+s2(3); 3.6+s2(4)],4);
A_new = vpa(x_new(1)*x_new(4) + x_new(2)*x_new(4)- x_new(1)*x_new(3),3);
I_new = vpa(x_new(1)*(x_new(4))^3 + x_new(2)*(x_new(4))^3 - x_new(1)*(x_new(3))^3 /12,3);
w_new = vpa(P_load + rho*A_new,3);
W_new = vpa(rho*L*A_new,3);
delta_new = vpa((w_new *(L^4))/(8*E*I_new ),3)
freq_new = vpa((3.52/(pi*2))*(sqrt((E*I_new*386.4)/(rho*A_new*L^4))),3)