%% Example 6.4
% Poisson Equation with non Homogeneous Boundary Conditions




%% Initialization
close all;
clear M N dx dy x y X Y Pji Pij Pjk fun Qij Qkj Qjk
clear k A b Sol

M = 32;
N = 32;
dx= 1/M;
dy= 1/N;
x = [0:dx:1];
y = [0:dy:1];
[X,Y] = meshgrid(x,y);

Pji = zeros(N+1,M+1);
Pij = zeros(M+1,N+1);
Pjk = zeros(M+1,N+1);


%% Define the RHS of modified eq
fun = @(x,y) 30*(x.^2-x) + 30*(y.^2-y) - 4*pi^2*(x-1).*sin(2*pi*y);
Qij = fun(X,Y);
Qkj = dst1(Qij');
Qjk = Qkj';


%% Solve using DST
for k = [1:M-1]
    %Solve the tridiagonal system for each k
    A = diag(ones(M-2,1),-1) + diag(ones(M-2,1),1) + ...
        diag((((dy/dx)^2*(2*cos(pi*k/M)-2))-2)*ones(M-1,1));
    A = sparse(A);
    b = dy^2 * Qjk(2:M,k+1);
    Pjk(2:M,k+1) = (A^-1)*b;
end;

%Inverse transform
Pji = dst1(Pjk')*M/2;
Pij = Pji';

%% Solution of the original equation
Sol = Pij - (X-1).*sin(2*pi*Y);


%% Generates the plots
surf(X,Y,Sol);
xlabel('x');ylabel('y');zlabel('phi');
title('Numerical solution of Poisson equation, ex 6.4');
