% This m-file answers the questions on plotting in Practical #3.
% Nick Hale, Asgeir Birkisson, Mohsin Javed

%% Problem 1

f=@(x) x.^10-0.1;
a=0; b=1;

% Define a tolerance:
tol = 1e-10;

% Initialize iteration counter:
nIter = 1;
fr=[];

% Loop and bisect:
while( abs(b-a) > tol && nIter < 100 )
    r = (a + b)/2;        
    fr(nIter)=f(r);
    
    % if root is in [a, r]
    if ( f(a)*f(r) <= 0 )
        % make r as the new b:
        b = r;
        % bisect again:
        r = (a + b)/2;
    end
    
    % if root is in [r, b]:
    if ( f(r)*f(b) <= 0 )
        % make r as the new a:
        a = r;
        % bisect again
        r = (a + b)/2;
    end
    nIter = nIter + 1;
end

if ( nIter >= 100 )
    r = NaN;
end

subplot(2,1,1), plot(1:nIter-1,fr,'x-')
xlabel('Iteration number'), ylabel('Residual')
% hard to see convergence after about 5 iterations, use log scale
subplot(2,1,2), semilogy(1:nIter-1,abs(fr),'x-')
xlabel('Iteration number'), ylabel('Residual')
            
%% Problem 2

%% 2a
close all
figure
[xx, yy] = meshgrid(-3:.1:3,-3:.1:3);
f1 = @(x, y) (x.^2 + 3*y.^2).*exp(-x.^2 - y.^2);
surf(xx,yy,f1(xx,yy))

%% An alterative is to use fsurf
fsurf(@(x,y)f1(x,y),[-3,3,-3,3])

%%  2b
figure
[xx, yy] = meshgrid(-3:.1:3,-4:.1:4);
f2 = @(x, y) -3*y./(x.^2 + y.^2 + 1);
surf(xx, yy, f2(xx, yy), 'EdgeAlpha', 0)

%% 2c
figure
[xx, yy] = meshgrid(-1:.05:1, -1:.01:1);
f3 = @(x, y) abs(x) + abs(y);
surf(xx, yy, f3(xx, yy), 'EdgeAlpha', 0.5)

%% 2d
%
% MATLAB does not do any plotting for entries which are equal to NaN -- 'not a
% number'. Thus, by replacing entries in the ZZ matrix outside of the disc with
% NaNs, they won't be plotted:
figure
[xx, yy] = meshgrid(-1:.05:1, -1:.01:1);
f4 = @(x,y) sin(x).*cos(y);
zz = f4(xx, yy);
zz(xx.^2 + yy.^2 > 1) = NaN;
surf(xx, yy, zz)
axis([-1 1 -1 1])

%%
% An alternative approach is to use a different coordinate system, going from
% polar to cartesian coordinates.
figure
h = 2*pi/100;
[rr, tt] = meshgrid(0:.05:1, 0:h:(2*pi+h));
[xx, yy] = pol2cart(tt, rr);
f4 = @(x,y) sin(x).*cos(y);
zz = f4(xx, yy);
surf(xx, yy, zz)
axis([-1 1 -1 1])

%% Problem 3
close all

figure
[xx, yy] = meshgrid(-3:.1:3, -3:.1:3);
f1 = @(x, y) (x.^2 + 3*y.^2).*exp(-x.^2 - y.^2);
mesh(xx, yy, f1(xx, yy), 'LineWidth', 1.1)
figure
fmesh(@(x,y)f1(x,y),[-3,3,-3,3], 'LineWidth', 1.1)

figure
contour(xx, yy, f1(xx, yy), 'LineWidth', 2)
figure
contourf(xx, yy, f1(xx, yy), 'LineWidth', 2)
figure
fcontour(@(x,y)f1(x,y),[-3,3,-3,3])