stupid matlab

beamy_function.m

function [ output_arg ] = beamy_function( x )
global r
output_arg = x^2*(1+(2*r*x)^2)^0.5;
end

comptrap.m

% Composite trapeziod rule
% ff - integrand function
%   a,b interval with b>a
% Useage
%  f=@(x) x; Comptrap(f,0,1,n)


function [sum]= Comptrap(f,a,b,n)
fprintf('Im in Comptrap\n')
sum=0;
h=(b-a)/n;
for v=1:n-1
x=a+h*v;
sum=sum+feval(f,x);
end
sum=h*(feval(f,a)+feval(f,b))/2+h*sum;
end

T3.m

%Calls comptrap, doubles n each time
n=1;
old=0;
a=0;
b=1;
i=0;
A=zeros(20,1);
B=zeros(20,1);

global r
        
for r=0:0.1:2
    
    i=i+1;
    
    for U=1:30
    
        fprintf('U = %d ',U)
        fprintf('Im here\n')

        f=@beamy_function; 
        ans1 = Comptrap(f,a,b,n);
        if abs(old-ans1)<0.005
            break
        end
        n=n*2;
        old=ans1;
    end
    
    dP=(10^3/10^4)*ans1;
    A(i)=dP;
    B(i)=r;
end

plot (A,B)