function OrtogonalitaLegendrovychPolynomu()

syms x
w = 1;

f = sin(x);
fn = 0;

f1 = figure;
f2 = figure;
xx = -1:0.1:1; plot(xx,eval(subs(f,xx)),'r')
N = 5;
for k=0:N
    if k==0
        P{1} = 1;
    elseif k==1
        P{2} = x;
    else
        P{k+1} = (2*(k-1)+1)/((k-1)+1)*x*P{k}-((k-1)/((k-1)+1))*P{k-1};
    end
    for j=0:k
        log = sprintf('(P_%d,P_%d)=%g',j,k,eval(int(w*P{j+1}*P{k+1},x,-1,1)));
        disp(log);
    end
    % Legendrovy polynomy
    figure(f1); hold on;
    yy = eval(subs(P{k+1},xx));
    plot(xx,yy);
    
    % aproximace sin(pi*x)
    figure(f2); hold on;
    alfa(k+1) = int(P{k+1}*f,x,-1,1)/int(P{k+1}*P{k+1},x,-1,1);
    fn = fn + alfa(k+1)*P{k+1};
    yy = eval(subs(fn,xx));
    plot(xx,yy);
    pause
end

