clear all
disp(' ')
n = input('Ge antal gridpunkter (obekanta), n = ');
nx=sqrt(n);
ny=nx;
N = nx*ny;
h = 1 / (nx+1);

%  Set up matrix A and vector f. Use MATLAB's sparse storage.

Adiag = 4*ones(N,1);
Asubs = -ones(N,1);
Asubs1 = -ones(N,1);
Asuper1 = -ones(N,1);

% we must account for the zeroes in the sub and super diagonal

for i=1:ny
  Asub1(i*nx) = 0; 
end

for i=0:ny-1
  Asuper1(i*nx+1)=0;
end

A = spdiags([Asubs,Asubs1,Adiag,Asuper1,Asubs],[-nx -1 0 1 nx],N,N);

A=A./h^2;



[x,y]=meshgrid(h:h:1-h,h:h:1-h);

f=zeros(nx,ny);
f=1000*exp(-1000*((x-0.5).^2+(y-0.5).^2));

F=reshape(f,nx*ny,1);

tic;
  solution=A\(F);
t1=toc;
disp(' ')
disp('RESULTAT')
disp(' ')
disp(['Uppmatt tid for att losa systemet, T = ' num2str(t1,'%12.6f') ' s'])
disp(' ')
disp(' ')

u=reshape(solution,nx,ny);

[XX,YY]=meshgrid(h:h:1-h,h:h:1-h);
figure, 
surf(XX,YY,u), shading interp
xlabel('x')
ylabel('y')
title('Temperaturfordelning over en platta')
colorbar

B=A;