format long
M_s = 0;
alpha_s = 0;
Theta_F = [1:1:15]
iMax = 5;
A1 = [0,-30,-120,-270,-480,-750,-1080,-1470,-1920,-2430,-3000,-3630,-4320,-5070,-5880]
B1 = [-196200, -196200, -196200, -196200, -196200, -196200, -196200, -196200, -196200, -196200, -196200, -196200, -196200, -196200, -196200]
f = @(x) A1.*cos(x) + B1.*sin(x) + M_s*(x + (Theta_F ./(180*pi)) + (alpha_s /(180*pi)));
fd = @(x) -A1.*sin(x) + B1.*cos(x) - M_s;
x0 = 0.0001;
x = zeros(1, length(A1));
x(:) = x0;
xSave = zeros(iMax, length(A1));
for i = 1:iMax
x = x - f(x)./fd(x);
xSave(i,:) = x;
i = i+1;
end
y = f(x)