%
% 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)

% Actually, it is not easy to build a matrix out of block matrices
% Or does anybody know how. (Compare Maple's band-command)

% Anyway, it is more efficient to form it diagonal by diagonal.

% 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=   % help ones
neardiag= ...
fardiag=  ...

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.

% help rem
% 1:N-1   % gives [1,2,...,N-1]
% use ==  or ~=  % Note: scalar extension and comparison componentwise

neardiag=...

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,0,1,n],N,N)
full(A)

% If you want to see a part of it type
% AF=full(A);
% AF(1:5,1:5)  % for instance
%Now, let's write it as a function:

%Save it in the file lapm.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Problem nr. 2
%%%%%%%%%%%%%%%
%
left=[20,40,60,80];top=[50,20,10,0];bot=zeros(1,4),right=left

m=5;n=m-1;

%First build B as a 4x4 matrix. B(:) makes it a vector (columnwise).
%Thus you need to transpose B before "raveling" it into a long vector.


A=lapm(m); % if you did'n write it as a function, cut and paste the commands
from above.

Tvek= ... Use \  (backslash)
T=reshape(Tvek,4,4)'  % reshape works columnwise, thus you need to transpose
T=flipud(T)           % Flip to fit the coordinate-axes
T=[fliplr(left)' T fliplr(right)']  % Vertical boundaries need flips as well

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));


% When you are done, check it!

for i=2:n; for j=2:n ;
  tt(i,j)= 4*T(i,j)-T(i-1,j)-T(i+1,j)-T(i,j-1)-T(i,j+1);
end; end
tt

% There's a function for this:

del2(T)