PROBLEM SET 1

% Exercise 1.1
% Demonstrates the relationship between Euler and trapz errors for multiple
% dx values
function errortest

dx_values = zeros(1,997);
euler_error = zeros(1,997);
trapz_error = zeros(1,997);

for i = 3:1000
    x = linspace(0,pi/2, i);
    dx = x(2) - x(1);
    dx_values(i-2) = dx;
    y = sin(x);
    euler_error(i-2) = 1-sum(dx * y(1:end-1));
    trapz_error(i-2) = 1-trapz(x,y);
end

plot(dx_values, euler_error, 'color','red'); hold on;
plot(dx_values, trapz_error, 'color', 'blue');
legend('Euler', 'trapz');
xlabel('dx'); ylabel('error');
title('Euler vs trapz error as a function of dx');



% Exercise 1.7
function gammatest

fprintf('factorial(4) = %g\n',factorial(4));
fprintf('gamma(5) = %g\n',gamma(5));
fprintf('factorial(.5) = ERROR!\n');%, factorial(0.5));
fprintf('gamma(1.5) = %g\n', gamma(1.5));

end


% Exercise 1.5
function kevin

x = 0:pi/100:2*pi;

figure; hold on;
%for i = 1:10
    y1 = expandcos(x,7);
    y2 = expandsin(x,7);
    Y1 = sin(x);
    Y2 = cos(x);
    plot(x,y1,'b');
    plot(x,y2,'r');
%end
plot(x,Y1,'g'); plot(x, Y2, 'g');
axis([0 2*pi -1 1]);
end

function val = expandsin(x, n)
    parity = 1;
    val = 0;
    for i = 0:n
        val = val + x.^(2*i+1)/factorial(2*i+1)*parity;
        parity = -parity;
    end
end

function val = expandcos(x, n)
    parity = 1;
    val = 0;
    for i = 0:n
        val = val + x.^(2*i)/factorial(2*i)*parity;
        parity = -parity;
    end
end

-----------------
PROBLEM SET 2

% Exercise 2.1
function z = atantest(y,x)
if(y > 0)
    z = atan2(y,x);
else
    z = atan2(y,x)+ 2*pi;
end
end

2.3/2.4 - can't figure out how to use ODE23


% Exercise 2.5
function einstein(filename)
load(filename);
hist(einstein,100);
end%

Exercise 2.6
Y = inline('sin(pi*x).^2');
quad(Y,0,1)

% Exercise 2.7
function twopointthreeseven
x = rand(1,1e4);
y = cheese(x);

%plot(x,y);
[N, X] = hist(y,25);
area = sum(N.*10e-4.*X);
fprintf('Area calculated from histogram: %g\n', area);


function y = cheese(x)
    y = sin(pi*x).^2;
end

end

% Exercise 2.8
function twopointeight

x = rand(1,1e4);
y = rand(1,1e4);
z = x.^2+y.^2;

num_inside = length(find(z<=1));
fprintf('area of circle approximated: %g\n',num_inside*4/1e4);

end