% definizione problema differenziale
funp = @(x,y) x-1;
funq = @(x,y) 1-y;
funr = @(x,y) x-y;
funs = @(x,y) y-x+x.*y.^2-x.^2.*y;
a = -3;
b = 3;
c = -2;
d = 2;
funfa = @(y) -3*y;
funfb = @(y) 3*y;
fungc = @(x) -2*x;
fungd = @(x) 2*x;
funsol = @(x,y) x.*y;

% parametri e discretizzazione
n = 40;
m = 30;
h = (b-a)/(n+1);
k = (d-c)/(m+1);
x = linspace(a,b,n+2)';
y = linspace(c,d,m+2)';
[XG YG] = meshgrid(x,y);
X = XG(2:m+1,2:n+1);
Y = YG(2:m+1,2:n+1);

% visualizzazione soluzione
figure(1)
sol = funsol(XG,YG);
surf(XG,YG,sol)

% calcolo coefficienti sulla griglia
AM = funr(X,Y) - 2/(h^2) - 2/(k^2);
BM = 1/(h^2) - funp(X,Y)/(2*h);
BHM = 1/(h^2) + funp(X,Y)/(2*h);
CM = 1/(k^2) - funq(X,Y)/(2*k);
CHM = 1/(k^2) + funq(X,Y)/(2*k);
SM = funs(X,Y);

% calcolo delle condizioni al contorno
fa = funfa(y(2:m+1));
fb = funfb(y(2:m+1));
gc = fungc(x(2:n+1));
gd = fungd(x(2:n+1));

% trasposizione per rispettare la notazione del libro
AM = AM';
BM = BM';
BHM = BHM';
CM = CM';
CHM = CHM';
SM = SM';

% inserimento degli zeri in b e bhat
b1v = BM(1,:);
BM(1,:) = 0;
bhnv = BHM(n,:);
BHM(n,:) = 0;

% vettorializzazione dei coefficienti
av = AM(:);
bv = BM(:);
bhv = BHM(:);
cv = CM(:);
chv = CHM(:);
sv = SM(:);

% preparazione delle diagonali (shift)
bv = [bv(2:end);0];
bhv = [0;bhv(1:end-1)];
cv = [cv(n+1:end);zeros(n,1)];
chv = [zeros(n,1);chv(1:end-n)];

% matrice dei coefficienti
DIAGS = [cv bv av bhv chv];
N = n*m;
A = spdiags( DIAGS, [-n -1 0 1 n], N, N);

figure(2)
spy(A)

% termine noto
f = -sv;
f(1:n) = f(1:n) - cv(1:n).*gc;
f(end-n+1:end) = f(end-n+1:end) - chv(end-n+1:end).*gd;
f(1:n:end) = f(1:n:end) - b1v'.*fa;
f(n:n:end) = f(n:n:end) - bhnv'.*fb;

% sistema lineare
u = A\f;
U = reshape(u,n,m);
fa = funfa(y);
fb = funfb(y);
U = [fa'; gc U gd; fb'];
U = U';

figure(3)
surf(x,y,U);
title('approssimazione')

figure(4)
surf(x,y,sol-U)
title('errore')


