(*Implementacja kwadratury Gaussa z możliwością przeskalowania przedziału całkowania do (-1,1)*)
(*aa i bb to granice żądanego przedziału*)
aa=-1;
bb=1;
c1=(bb-aa)/2;
c2=(bb+aa)/2;
(*r[x] to funkcja wagowa*)
wiel[n_,x_]:=LegendreP[n,x];
r[x_]:=1;
a=-1;
b=1;

n=7;

absx=SetPrecision[NSolve[wiel[n,x]==0,x],50]

A=Table[wiel[w,x]/.absx[[p]],{w,0,n-1},{p,1,n}];

Norma=Integrate[wiel[0,x]^2*r[x],{x,a,b}]

V=ConstantArray[0,n];

V[[1]]=Norma/wiel[0,x];

w=LinearSolve[A,V]

f[x_]:=x^2;
(*UWAGA! "Ruszamy" tylko zmienną w funkcji f. Funkcja wagowa r[x] musi być "rozciągnięta" wraz z przedziałem całkowania*)

c1*N[Integrate[r[x]*f[c1*x+c2],{x,a,b}],50]

c1*Sum[w[[i]]*f[c1*x+c2]/.absx[[i]],{i,1,n}]

(*kopia punktów próbkowania funkcji, będzie dalej potrzebna*)

abst=absx/.x->t

(*POCZĄTEK KODU ROZWIĄZUJĄCEGO RÓWNANIA CAŁKOWE FREDHOLMA DRUGIEGO RODZAJU*)
(* f(x)=g(x)+Alfa*Integrate[K[x,t]*f[t]dt,{t,a,b}]*)
(*K - jądro*)

K[x_,t_]:=Sin[Pi*x*t/2];

Alfa=1/2;

(*Metoda Nystroma*)

g[x_]:=x^3;

Gvec=Table[g[x]/.absx[[i]],{i,1,n}]

MatW=IdentityMatrix[n]*w;

MatK=Table[K[x,t]/.absx[[i]]/.abst[[j]],{i,1,n},{j,1,n}];

wynik=LinearSolve[(IdentityMatrix[n]-Alfa*MatK.MatW),Gvec]

Funkcja=g[x]+Alfa*Sum[w[[i]]*(K[x,t]/.abst[[i]])*wynik[[i]],{i,1,n}]

Plot[Funkcja,{x,-1,1}]

(*Metoda kolejnych przybliżeń*)

Funkcja=1

Do[Całka=Sum[w[[i]]*(K[x,t]/.abst[[i]])*(Funkcja/.absx[[i]]),{i,1,n}];
Funkcja=g[x]+Alfa*Całka;
,{j,1,10}]
Funkcja

Plot[Funkcja,{x,-1,1}]