% Solutions # 1
%
% Problem 1
% 
n=3; B=4*diag(ones(1,n))-diag(ones(1,n-1),1)-diag(ones(1,n-1),-1)
O=zeros(n,n)
mI=-eye(n,n)

% Let m be the nr. of subintervals
% n=m-1   % nr. of interior grid point in each axis
% N=N^2   % Nr. of inner nodes on the 2-dim. grid = nr. of unknowns

m=4;n=m-1;N=n^2
maind=4*ones(N,1)
neardiag=-ones(N-1,1);
fardiag=-ones(N-n,1)

diag(maind)
diag(neardiag,1)
diag(neardiag,-1)
diag(fardiag,n)
diag(fardiag,-n)

diag(maind)+...
diag(neardiag,1) +...
diag(neardiag,-1) + ...
diag(fardiag,n) +...
diag(fardiag,-n)

% It's good, except you need to have 0 at each nth position of neardiag.

neardiag=rem(1:N-1,n)~=0

A=...
diag(maind)+...
diag(neardiag,1) +...
diag(neardiag,-1) + ...
diag(fardiag,n) +...
diag(fardiag,-n)


%%%%%%%%%%%%%%%%%%
% Sparse matrices:
help sparse
help spdiags

% The easiest way to convert the matrix in sparse storage is
A=sparse(A)

% This is not an economic way to build a sparse matrix though. 
% (Often the full matrix doesn't fit into memory)
%
% spdiags offers a way to build banded sparse matices storing just the
% diagonals.
%
m=4;n=m-1;N=n^2
maind=4*ones(N,1)
fardiag=-ones(N,1)
neardiag=rem(1:N,n)~=0
[fardiag,neardiag',maind,neardiag',fardiag]

% This looks good but isn't quite. There's a little trick about spdiags
% not explained in help spdiags. It uses the data 
%   starting at the top on the LHS (below) the maindiag but
%   starting at the bottom on the RHS (above) the maindiag.
% (This doesn't show on constant diagonals.)

% Thus we need to append 1 at the head and drop one extra element at the
% tail of the upper neardiag.


uneardiag=[1,neardiag];uneardiag=uneardiag(1:N);
[fardiag,neardiag',maind,uneardiag',fardiag]
A=spdiags([fardiag,neardiag',maind,uneardiag',fardiag],[-(n+1),-1,0,1,n+1],N,N)
full(A)

%Now, let's write it as a function:

%Save it in the file lapm.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function A=lapm(m)
% Matrix for solving Laplace (Poisson) equation with 5 point. discr.
% Sparse structure.
% Arguments: m -- nr. of subintervals
%
n=m-1;  % 
N=n^2; % nr. of inner nodes = nr. of unknowns
e=ones(N,1);
d=~rem(1:N,n)==0;du=[1,d];du=du(1:N);
A=spdiags([e,d',-4*e,du',e],[-(n+1),-1,0,1,n+1],N,N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Problem nr. 2
%%%%%%%%%%%%%%%
%
left=[20,40,60,80];top=[50,20,10,0];bot=zeros(1,4),right=left

m=5;n=m-1;
Z=zeros(1,n-2);
B=[bot;[left(2:n-1)' zeros(n-2,n-2) right(2:n-1)'];top]
B(1,1)=B(1,1)+left(1);B(1,n)=B(1,n)+right(1);
B(n,1)=B(n,1)+left(n);B(n,n)=B(n,n)+right(n);
B=B';B=B(:)


A=-lapm(m);

Tvek=A\B
T=reshape(Tvek,4,4)'
T=flipud(T)

T=[fliplr(left)' T fliplr(right)']

NW=0.5*(left(n)+top(1));NE=0.5*(top(n)+right(n));
SE=0.5*(bot(n)+right(1));SW=0.5*(bot(1)+left(1));

T=[NW top NE;
   T;
   SW bot SE]

del2(T)

