I don't know if this is what you want. Your code is very badly structured and hard to follow.
Shouldn't
x=[2 4 6 8 10 12 14 16 18]; %x diameter in (cm)
instead of
x=[2 4 6 8 10 12 14 16 8]; %x diameter in (cm)
at the two places in your code ?
clear, close
%%%%NEWTON DIVIDED DIFFERENCE%%%
%all numbers in (cm)
l=50; %length
x=[2 4 6 8 10 12 14 16 8]; %x diameter in (cm)
fx=[5 7 8.4 9.1 9.1 8.3 8.1 7.8 6.6];%y diameter in (cm)
r=7; %radius of final bar
%INPUT Newtons divided difference
k=length(x);
B=zeros(k,k);
B(:,1)=fx;
for j=1:k-1
for i=1:k-j
B(i,j+1)=(B(i+1,j)-B(i,j))/(x(i+j)-x(i));
end
end
syms X
f=0;
m=1;
f=f+B(1,1)*m;
for i=1:k-1
m=m*(X-x(i));
f=f+B(1,i+1)*m;
end
p=3.3 ;
A=double(subs(f,X,p));
display(A)
syms U
f=0;
m=1;
f=f+B(1,1)*m;
for i=1:k-1
m=m*(U-x(i));
f=f+B(1,i+1)*m;
end
p=2 ;
A=double(subs(f,U,x));
display(A)
syms X;
f=B(1);
for i=1:k-1
f=f+B(i+1)*X^i; %f=required polynomial
end
disp(f);
syms X x
f=subs(f,X,x);
%Trapezodial Rule
%INPUT SECTION
f = matlabFunction(f);
f = @(x)f(x)-0.5*pi*r^2;
%f=@(x)f-(0.5*(pi*r^2)); %function
x=[2 4 6 8 10 12 14 16 8];
a=x(:,1); %lower limit
b=x(:,k); %upper limit
tol=1/100; %tolerance= percentage given/100
I0=0;
RE=100;
%CALCULATIONS
n=0;
while RE>tol
n=n+1;
h=(b-a)/n;
x_vec=a:h:b;
for i=1:length(x_vec)
fxx(i)=f(x_vec(i));
end
I=(h/2)*(fxx(1)+2*sum(fxx(2:n))+fxx(n+1));
RE=abs((I-I0)/I)*100;
fprintf('n1=%d \t I=%0.4f\t RE=%0.2f \n',[n;I;RE])
I0=I;
end
