np = 300;
rngseed = 1;
xmin = 0;
xmax = 5;
[x, ynoisy, yclean] = noisydata(@fun, xmin, xmax, np, rngseed);

plot(x, yclean, 'k--');

hold on

plot(x, ynoisy, 'b.')

A = [phi1(x') phi2(x') phi3(x')];
c = myqrfit(A, ynoisy)
yfit = c(1)*phi1(x) + c(2)*phi2(x) + c(3)*phi3(x);

plot(x, yfit, 'r-')

grid on
legend('Original function', 'Noisy data', 'LSQ fit')
xlabel('x')
ylabel('y')
title('Least squares fit')
hold off

function f = phi1(x)
    f = x;
end

function f = phi2(x)
    f = cos(x);
end

function f = phi3(x)
    f = exp(-x);
end
    
function y = fun(x)
    a = 1.;
    b = 1.;
    c = 1.;
    y = a*phi1(x) + b*phi2(x) + c*phi3(x);
end

